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Abstract. Efficient implementations of electronic structure methods are essential 
for first-principles modeling of molecules and solids. We here present a particularly 
efficient common framework for methods beyond semilocal density-functional theory, 
including Hartree-Fock (HF), hybrid density functionals, random- phase approximation 
(RPA), second-order M0ller-Plesset perturbation theory (MP2), and the GW method. 
This computational framework allows us to use compact and accurate numeric atom- 
centered orbitals (popular in many implementations of semilocal density-functional 
theory) as basis functions. The essence of our framework is to employ the "resolution of 
identity (RI)" technique to facilitate the treatment of both the two-electron Coulomb 
repulsion integrals (required in all these approaches) as well as the linear density- 
response function (required for RPA and GW). This is possible because these 
quantities can be expressed in terms of products of single-particle basis functions, which 
can in turn be expanded in a set of auxiliary basis functions (ABFs) . The construction 
of ABFs lies at the heart of the RI technique, and here we propose a simple prescription 
for constructing the ABFs which can be applied regardless of whether the underlying 
radial functions have a specific analytical shape (e.g., Gaussian) or are numerically 
tabulated. We demonstrate the accuracy of our RI implementation for Gaussian and 
NAO basis functions, as well as the convergence behavior of our NAO basis sets for 
the above-mentioned methods. Benchmark results are presented for the ionization 
energies of 50 selected atoms and molecules from the G2 ion test set as obtained with 
GW and MP2 self-energy methods, and the G2-I atomization energies as well as the 
S22 molecular interaction energies as obtained with the RPA method. 
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1. Introduction 



Accurate quantum-mechanical predictions of the properties of molecules and materials 
(solids, surfaces, nano-structures, etc.) from first principles play an essential role 
in chemistry and condensed-matter research today. Of particular importance are 
computational approximations to the many-body Schrodinger or Dirac equations that 
are tractable and yet retain quantitatively reliable atomic-scale information about the 
system — if not for all possible materials and properties, then at least for a relevant 
subset. 

Density-functional theory (DFT) [TJ |2] is one such successful avenue. It maps 
the interacting many-body problem onto an effective single-particle one where the 
many-body complexity is hidden in the unknown exchange-correlation (XC) term, 
which has to be approximated in practice. Existing approximations of the XC term 
roughly fall into a hierarchial scheme [3]. Its local-density (LDA) [2 J and generalized 
gradient approximations (GGAs) [U EJ EJ [7] are now well recognized workhorses with a 
broad application range in computational molecular and materials science. However, 
several qualitative failures are well known: To name but a few, certain adsorbate 
geometries [SJ, /-electron systems 0, EH EH U21 E3], or van der Waals interactions 
[HI EH EH El EH El] are not described correctly at this level of theory. Thus, there is 
much ongoing work to extend the reach of density-functional theory, e.g., meta-GGAs 
[201 HH [22], formalisms to include van der Waals interactions [2^ [2^ [25] I26[ |27[ 128]. 
hybrid functionals (29j [30l EH E21 EHl [34], or approaches based on the random-phase 
approximation (RPA) [35l ESI EZl [381 E21 SOI SH S2] that deal with the non-local 
correlations in a more systematical and non-empirical way. 

Another avenue are the approaches of quantum chemistry, that start with Hartree- 
Fock theory [431 EI] • These approaches offer a systematically convergable hierarchy of 
methods by construction, said to reach "gold standard" accuracy for many molecular 
systems at the level of coupled-cluster theory [451 US E] [often, taken to include singles, 
doubles, and perturbative triples, CCSD(T) [48]]. CCSD(T) theory is significantly more 
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accurate than DFT-LDA/GGA for many molecular systems but also significantly more 
costly (formally scaling as 0(N 7 ) with system size). It has its own shortcomings as 
well. For instance, systematic, material-specific failures can occur in cases where the 
underlying Hartree-Fock solution itself is not a good reference to start with (for example, 
many open-shell systems), and a multireference extension of the approach [19] becomes 
necessary. 

A third avenue for electronic structure calculations is the quantum Monte Carlo 
(QMC) method, in particular the diffusion QMC method [50], EI]. This is a stochastic 
approach that deals with the many-body wavefunction directly. The diffusion QMC 
method can often deliver high accuracies, and provide data and insights for problems 
which are difficult for other approaches. Its widespread use, however, is also impeded by 
the rather high computational costs. Moreover, the fixed-node approximation and the 
underlying pseudopotential approximation are known issues which limit the practical 
accuracy of QMC. Regarding the computational cost, the QMC methods, the algorithm 
of which is intrinsically parallel, are in a better position to benefit from the development 
of petaflop supercomputers [52] . 

With the successes, but also the failures or shortcomings of the aforementioned 
avenues, much attention is currently devoted to the construction of further, systematic 
and generally applicable methods or theoretical frameworks that can offer better 
accuracies than conventional DFT, but have lower numerical costs and are free of the 
limitations of CCSD(T) and QMC. Among the various possible pathways, many-body 
perturbation theory (MBPT) based on an efficiently attainable and trustful electronic 
reference state offers such an avenue. In particular, approaches based on the RPA, 
which bridge the DFT and MBPT worlds [35], [36] El] HE], have recently enjoyed 
considerable attention for ground-state total-energy calculations. For electron addition 
and removal energies, a self-energy based approach that is consistent with the RPA 
total-energy treatment is Hedin's GW approximation [M]- GW is especially popular 
in the solid state community [55] EE] EH EH] and has become the method of choice 
for the calculation of quasiparticle band structures as measured in direct and inverse 
photoemission [59] [60], EI] . 

Although RPA and GW are receiving much attention in the community today, 
the systematic investigation of diagrammatic perturbation theory from first principles 
for real materials is only just beginning. Its full promise lies in the fact that it 
is intermediate in cost between DFT and coupled-cluster theories, and applicable in 
practice to molecular and condensed materials alike - including open-shell systems and 
metals. 

Besides the more generally applicable RPA and GW approaches, another correlation 
method that is widely used in computational chemistry is second-order M0ller-Plesset 
theory (MP2) [621 HI], which belongs to the category of Hartree-Fock based quantum 
chemistry approaches mentioned above. MP2 does not reach the CCSD(T) accuracy, 
but its more favorable computational scaling makes it applicable to larger system 
sizes. In analogy to the GW self-energy, a MP2 self-energy approach [63] 03] that 
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is compatible with the MP2 total energy is possible as well. As will be demonstrated 
more clearly in the next section, RPA, GW, and MP2 (both total and self-energy) are 
related approaches both diagramatically and numerically. The development of numerical 
frameworks that enables the implementation of all these approaches on an equal footing 
with high numerical efficiency and accuracy is thus highly desirable. 

In the present work, we describe the underpinnings of such a unified numerical 
framework that is promising to boost the efficiency for all the above-mentioned methods, 
by allowing for their implementations with compact and efficient NAO basis sets. Our 
specific implementation is based on the "Fritz Haber Institute ab initio molecular 
simulations" (FHI-aims) [61] program package. While we make reference to FHI-aims 
basis sets throughout much of this work, the numerical foundation presented here is 
general: applicable to any other type of atom-centered basis set. We note that many 
production-quality implementations of hybrid functionals, Hartree-Fock, and MBPT 
are based on analytical basis functions such as Gaussian-type orbitals (GTOs) or plane 
waves, and typically rely on pseudopotential-type approaches. In contrast, we here aim 
for an all-electron, full-potential treatment with NAO basis sets that does not sacrifice 
accuracy compared to the alternatives. 

For DFT-LDA/GGA, NAOs are well established and can be found in several 
implementations [65j ESI EJJ EHl EH EHJ E0]- This is however not the case for HF and the 
MBPT approaches we are going to address in this paper. Specifically we will present in 
the following 

i) an atom-centered resolution of identity (RI) framework analogous to what is 
pursued in quantum chemical methods [7H E21 E31 EH EH EH E7J EH]- This 
framework allows us to reduce all four-center two-electron Coulomb integrals to 
precomputed three- and two-center integrals. Our scheme differs from the quantum 
chemistry approach [TTJ EH] in the auxiliary basis set construction, which is essential 
for retaining the flexibility to work with any atom-centered basis function shape, 
rather than being restricted to analytical shapes only. 

ii) an assessment of the accuracy of the NAO basis sets used for normal LDA/GGA 
calculations in FHI-aims, and intended to be transferable regardless of the specific 
underlying materials or functionals, for Hartree-Fock, MP2, hybrid functionals, 
RPA, and GW. 

Reference to relevant work by other groups in electronic-structure theory is made 
throughout this work. 

The present paper demonstrates our approach for molecular systems (non-periodic), 
and makes extensive use of established GTO basis sets for comparison and reference 
purposes. In addition, we provide benchmark GW vertical (geometry of the ionized 
molecule not relaxed) ionization energies (IEs) for a subset of the G2 ion test set [79], 
and benchmark binding energies (RPA) for the G2-I and the S22 molecular test set ]80j . 
We restrict ourselves to algorithms that have standard scaling with system size [0(N 4 ) 
for HF, 0(N 5 ) for MP2, etc]. Regarding total energy differences we restrict ourselves to 
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a discussion of counterpoise-corrected [ST] results when assessing the accuracy of MBPT 
methods that utilize the full (also unoccupied) spectrum. 

In the following, we first recapitulate the HF, MP2, RPA, and GW methods 
and highlight the structural similarity and difference of the three correlation methods 
(section [2]) . We then introduce the basics of RI and the RI formulation of the above 
methods (section^. Our own RI prescription and its accuracy is the subject of section[U 
Section [5] demonstrates the overall accuracy of our approach with NAO basis sets for 
a variety of test systems for HF, MP2, hybrid density functionals, RPA, and GW. 
The benchmark results for a subset of the G2 and S22 molecular test sets using our 
approaches are presented in section El Finally we conclude our paper in section [7J 



2. Theoretical framework: HF, hybrid density functionals, MP2, RPA, and 
GW 



2.1. Many-electron Hamiltonian and many-body perturbation theory 

Hartree-Fock (HF) theory, hybrid density functionals, and MBPT (MP2, RPA, and 
GW) are all approximate ways to solve the interacting many-electron Hamiltonian 
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where N e is the number of electrons that interact via the Coulomb interaction vf* 



I r, — rj\, and vf* = u ex t( r i) is a local, multiplicative external potential, usually due 
to the nuclei. Hartree atomic units are used throughout this paper. The numerical 
cost for an exact solution of the Hamiltonian (CQ) scales exponentially with system size 
(number of electrons). The systems for which such a solution is possible are thus heavily 
restricted in size. In general, accurate approximations are needed. The most common 
approximations first resort to the solution of a mean-field, non-interacting Hamiltonian 
H° that yields an approximate ground-state wave function |$o) : 
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The "(0)" in Eq ^ implies the fact that this is the ground state energy of the mean- field 
Hamiltonian. A suitable H° should be solvable with relative ease. |$q) is a single Slater 
determinant formed from the lowest N e single-particle spin-orbitals determined by 



M>na) 



(3) 



where a denotes the spin index and h is the effective single-particle Hamiltonian noted 
in the bracket of (|2l) . The form of ([2J13]) is, of course, precisely that of Kohn-Sham(KS) 
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DFT (with a local mean- field potential fj MF ), or of HF theory and hybrid functional 
(with a non-local mean- field potential f 4 MF ). 

The purpose of starting with (jTJOJ) is to establish our notation for the following 
sections and to distinguish between 

• H (the many-electron Hamiltonian) , 

• the mean- field Hamiltonian H°, the solutions of which are many-electron wave 
functions given by single Slater determinants, and define an excited-state spectrum 
of their own (obviously not the same as that of H), 

• the effective single-particle Hamiltonian h°, which generates the single-particle 
orbitals ipna and orbital energies e na . 

In MBPT, one starts from h° and the associated eigenenergies and eigenfunctions to 
systematically approximate the properties of H, e.g., its true ground-state energy E in 
MP2 or RPA or its singple-particle excitations in GW. The interacting many-electron 
Hamiltonian H is partitioned into a mean-field Hamiltonian H° as given by (J2J) and an 
interaction term H' , 

H = H° + H' 

i<j ' 1 J '' i=l 

In the remainder of this section we collect the basic formulae that define the mean- 
field Hamiltonians, perturbation theory for ground state properties (MP2, RPA), and 
perturbation theory for excited states (electron addition and removal energies, through 
either GW or MP2 self-energies). From a numerical point of view, the underpinning 
of all these methods is the same: an efficient, accurate basis set prescription, and an 
efficient expansion of the two-electron Coulomb operator, which is the primary focus of 
this paper. 



2.2. Mean-field Hamiltonians of HF or DFT 



In HF theory the ground-state wave function of the Hamiltonian in ([T]) is approximated 
by a single Slater determinant |$o) an d -Eg '' is obtained by a variational optimization, 
leading to 

(r|/!Vn) 



1. 

2 



-^V 2 + v ext (r) + v h (r) 



Vw(r) + / dr'S^ry)^ 1 ") = e na ip na (r) . (5) 



/ here denotes the HF single-particle Hamiltonian, and v h (r) is the Hartree potential 

n(r' x 



with the electron density 



v h (r) 



n(r) 



r — r' 



T dr' 
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and E* is the non-local, exact-exchange potential 

g(r.O = -E *-|/y • < 8 ' 

Equations ^SJ)— (JSJ) form a set of non-linear equations that have to be solved iteratively. 
v h (r) and £*(r, r') together yield the HF potential v HF , a special case of the mean- 
field potential in (j2J) and (HI). The HF wavefunction |$o) is given by the Slater 
determinant formed by the N e spin-orbitals ip nu with lowest energies e na . 
At self-consistency, the HF total energy is 

occ 

E RF = (%\H° + H'\%)=J2 <W — E h — E x 

na 

where the Hartree energy E h and exact-exchange energy E x are given respectively by 
E h = — I drn(r)uh(r) 



EX =2^JJ ^ r V; CT (r)S^r,r')Vn CT (r') 

na 

= --f^ff irir ^^( r )^( r )^( r/ )^( r . (9) 



In general, the single-particle spin-orbitals ip na (r) are expanded in terms of a set of basis 
functions {y?j(r)} 

Vw(r)=^4^(r). (10) 

i 

where c % na are the expansion coefficients. In terms of these basis functions, the HF 
effective potential can be expressed in a matrix form 

VZ = jf dvdr'^v) [v h (v)5(r - r') + S?(r, r')] V] (v') 

= <. a + S^, (11) 
where in particular the exact-exchange matrix is given by 

A;/ 

In (|T2|) -Dfci j(T is the density matrix 



occ 



n 



and {ij\kl) is the short-hand notation of quantum chemistry for 4-center 2-electron 
integrals 

mi) =11 ^ (r) ^ r) ^' M(r W . (i4) 



Very similar equations arise in KS-DFT with a local t^ MF , or in the generalized KS 
scheme [82J with a fraction of £*(r, r') in the potential. In principle, the exact KS-DFT 
would yield the exact many-electron ground-state energy E and ground-state density 
n (r). In practice, the XC energy functional and potential have to be approximated. 
The effective single-particle orbitals from either HF or from approximate KS-DFT are 
convenient starting points for MBPT. 



2.3. Perturbation theory for the many- electron ground-state energy: MP2 



Assuming that H' is benign and can be treated as a perturbation on top of H°, 
the ground-state energy for the interacting system can be obtained using Rayleigh- 
Schrodinger perturbation theory (RSPT), or more precisely Brueckner-Goldstone 
perturbation theory [831 El]- According to the Goldstone theorem [M], in a 
diagrammatic expansion of the ground-state total energy, only the "linked" diagrams 
need to be taken into account. And this guarantees the the size-extensivity of the 
theory, i.e., the total energy scales correctly with the system size. M0ller-Plesset (MP) 
perturbation theory is a special case of RSPT [62], where the reference Hamiltonian H° 
the HF Hamiltonian iffjF — Yli fi- Terminating the expansion at second order gives the 
MP2 theory, in which the (second-order) correlation energy is given by 



E, 



(2) 



k>0 



[®k\H'\<5> )\ : 



E 



(o) 
o 



E. 



(o) 



(15) 



Here are the Slater determinants representing the excited states of H 







H HF , 

and E^' are the corresponding excited-state energies. H' is given by (HI) with 
v MF = f HF . Among all possible excited-state configurations only double excitations 
contributes in (|T5l) . This is because singly-excited &k do not couple to the ground- 
state $o (Brillouin's theorem [44J for the HF reference), whereas even higher-excited 
configerations (triples, quadruples, etc.) do not contribute due to the two-particle nature 
of the operator H' . As such, equation ( jl~5]) can be expressed in terms of single-particle 
spin-orbitals, 



occ unocc 



ab 



[am, a\bn, a') — (bm, a\an, a')S a 



— £ba' 



(16) 



where (ma,a\nb,a') are 2-electron Coulomb repulsion integrals for molecular orbitals 

^>)^ CT (r)^(r')^(r') 



[ma, o~\nb, a') 



alralv 



r — r 



(17) 



The two terms in (1161) correspond to the 2nd-order Coulomb (direct) and 2nd-order 
exchange energy, respectively. 
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2-4- Perturbation theory for the many- electron ground-state energy: RPA 

MP2 corresponds to the 2nd-order term in a perturbation theory where the perturbation 
expansion is essentially based on the bare Coulomb operator (with HF effective potential 
subtracted). As such, the MP2 correlation energy diverges for the homogeneous electron 
gas and metals with zero direct energy gap. To overcome this problem in the framework 
of MBPT it is essential to sum up the diverging terms in the perturbation series to 
infinite order. One such example, which has gained considerable popularity recently 
[381 EH E6l EZl EH [891 Eol ISU E21 [391 [931 EH ESI [961 [9TI [98] , is the RPA [351 HE Ell E31 US], 
that in the context of MBPT corresponds to an infinite summation of "ring" diagrams. 
The choice of the non-interacting reference Hamiltonian H° can be, e.g., HF or DFT 
with any desired XC functional. 

Apart from the diagrammatic representation, RPA can also be formulated in other 
ways, e.g., as the simplest time-dependent Hartree approximation in the context of the 
adiabatic-connection fluctuation-dissipation theorem [99l [53J [16], or as a subclass of 
terms in coupled-cluster theory with double excitations [HZ]. In the context of DFT 
[53J, RPA calculations can be performed self-consistently by means of the optimized 
effective potential approach pm IM1 [T02| fTOa]. 

In close-packed notation, the RPA correlation energy (cRPA) can be expressed 
in terms of the RPA dielectric function e, or alternatively the non-interacting density 
response function x° on the imaginary frequency axis 

E RPA = ±_ / duTT [ ln + (1 _ 

i r°° 

= — J dcuTr [In (1 - X °Mv) + xVH 

1 POO °° 1 

= -o- / ^£-Tr[(xV>)1 • ( 18 ) 

171 J n =2 U 

The real-space (Adler- Wiser |104[ 1105] ) representation of x° reads 
X°(r, r',iu) = (r|x°(iw)|r / ) 



occ unocc 



Yl ^* ma ^ aa ^* aa ^ ma ^ I cc (19) 



~aa i '-ma 

a m a 



where cc. denotes "complex conjugate", and ^(r) and e n are single-particle orbitals 
and orbital energies as implied by (JHJ). The RPA dielectric function e is linked to x° 
through 

e(r, r, ico) = 5{r - r') - J dr"v(r, r") X °(r", r', ico) . (20) 
Using (I17jl and ( Tl9l) . the lowest-order term in ( IT8l) can be expressed as 

„ „ occ unocc / I s 

Tr \ X °(iu)v\ = / / drdv' X °(v, r', iu)v(r', r) = V V V y ma ^\ am ^> + cx . . (2 i) 
J J ^ ^ iu - e aa + e ma 



a m a 
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Higher-order terms in (1181) follow analogously. There is thus a straightforward path to 
compute x° ; an d hence the RPA correlation energy, once the selected non- interacting 
reference state is solved. In practice the RPA correlation energy is always combined with 
the exact-exchange energy (EX) in ([9}, henceforth denoted (EX+cRPA), but evaluated 
with the same single-particle orbitals as used in Ef FA . The choice of input orbitals will 
in the following be marked by (EX+cRPA) @MF, where MF specifies the mean- field 
approach used to compute the single-particle orbitals. The application of EX+cRPA to 
various realistic systems, as well as the development of schemes that go beyond simple 
EX+cRPA, is an active field [381 ESI ESI [1061 DUO E3 EH ED EE M> M> M> E31 HO 
[951 1961 [971 [981 11091 HO] . For instance, when the so-called single-excitation and second- 
order screened exchange contributions [1101 [90] are added to EX+cRPA, the resulting 
accuracy is impressive [40], [41] . 



2.5. Perturbation theory for electron addition or removal energies: GW or MP2 

Ground-state energies aside, one is often interested in the properties of electronically 
excited states. Part of this information is in principle accessible by taking the difference 
of the total energies of iV-electron system and N ± 1-electron systems using approaches 
discussed above. In practice, this approach mainly works for computing core-level 
excitations and/or the first ionization energy and electron affinity for (small) finite 
systems. In contrast, Green function techniques are more convenient and powerful for 
dealing with electronic excitations in general. The basic theory of Green functions is well 
documented in textbooks [111[ 1112] . Here we collect the contextual equations, based on 
which practical approximations can be introduced. 

The single-particle Green function of an interacting many-electron system is defined 

as 



G(v,t;r',t') = -i(N\T4>(r,t)ft(r',t')\N) 



(22) 



where \N) = |\l/o( r i; • • • > t n)) denotes the interacting ground-state wave function of 
an iV-electron system (solution to (CQ)). ip(r,t) and ^(r', t') are field operators in the 
Heisenberg picture that annihilate and create an electron at space-time point (r, t) and 
(r',t'), respectively. T is the time-ordering operator. The Green function G(r,t;r',t') 
measures the probability amplitude of a hole created at (r, t) propagating to (r', t') for 
t < t', or an electron added at (r', t') propagating to (r, t) for t > t' . The poles of 
its Fourier transform, G(r, r', u), correspond to the single-particle excitation energies as 
measured for example in direct and inverse photoemission experiments. 

For Hamiltonians with a time-independent external potential, as considered in 
this paper, the Green function depends only on the difference between t and t', 
G(r, t; r', t') = G(r, r'; t — t'). Its Fourier transform gives the frequency-dependent Green 
function G(r, r'; oj) that satisfies the Dyson equation, 



G(r, r'; to) — J <ir"£(r, r", co)G(r", r'; tu) = 5(r 



')■ (23) 
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Here v h (r) is the electrostatic Hartree potential defined in (E]) and E(r, r", oj) is the 
dynamical, non-local, complex self-energy that contains all the many-body XC effects. 

G(r,r';u>) and E(r, r';oj) for many-body interacting systems can in principle be 
obtained using diagramatic Feynman-Dyson perturbation theory. The perturbation 
series is built on a non-interacting Green function G°(r,r',u) that corresponds to the 
non-interacting single-particle Hamiltonian h°. With single-particle orbitals ip na {r) and 
orbital energies e na of h°, one has 

G°(r,r» = £ <MrK '< r,) p (24) 

^ oj - e na - irj sgn(e F - e na ) 

where e-p is the Fermi energy, and 77 a positive infinitesimal. For a given G°, 
corrections to the single-particle excitation energies can be computed from approximate 
perturbative expansions of E(r, r'; oj). Examples are the GW method and the 2nd-order 
approximations discussed below. 

In the GW approximation proposed by Hedin [51] , the self-energy assumes the form 

Ef^(r, r', oj) = ^ J dojG a {v, r', oj + oj')W{v, r', oj')e^ . (25) 

Here W is the screened Coulomb potential at the RPA level 

W(r,r',co) = j ' dr"e- 1 {r,T , \uj)v{v'\Y) , (26) 

with the dynamical dielectric function e as defined in (Tl9|) and (!20l) . but on the real 
frequency axis with ioj replaced by oj + irj (i] — > + ). 

In practice, one-shot perturbative GW calculations (often referred to as G ^ ) 
based on a fixed, non-interacting reference state (DFT with popular functionals or HF) 
are often performed. With x° computed from the non-interacting Green function in 
(|24|) . the ensuing W can be expanded in powers of x° v 

W = v + vx°v + vx°vx°v H . (27) 

The GW approximation can thus be regarded as an infinite series in the bare Coulomb 
potential v , or alternatively, as is obvious from ( l25i) . as first-order perturbation in terms 
of the screened Coulomb potential W. 

Once the G°W° self-energy is obtained from ( 1251) . the corrections to the single- 
particle orbital energies are given by 

eZ W ° = e n « + (^f W ° ( e G n : w °) ~ v^ n „) , (28) 



where the XC part v xc of the reference mean-field potential has to be subtracted. 
1281 approximates the quasiparticle wavefunctions with the single-particle orbitals of 
the reference state. This is often justified, but may break down in certain cases 

[munninainEiiirziinB]. 
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Rewriting the diagonal elements of the G°W° self-energy on the imaginary 
frequency axis in terms of four-center Coulomb integrals, we obtain 

1 f°° 

= -— du'G° m(7 (iu + iuj')(nm,a\W (iu')\mn,a) (29) 

where G° ma (iuj) = l/(iu + e F — e mcr ), and 

(nm,a\W°(iuj)\mn,a) = J J drdr / ^(r)V' mff (r)W(r,r / ,ia;)^(r , )^(r / ) . (30) 

This expression can be analytically continued to the real-frequency axis using Pade 
approximation or a two-pole model [119] . 

The G°W° approach, as a highly popular choice for quasiparticle excitation 
calculations, in particular for solids [59 | I12Q|, [60] , [6T| , I121j . is akin to the RPA approach for 
computing the ground-state correlation energy. Similarly, one can also introduce a 2nd- 
order perturbation theory for the self-energy [Sj consistent with the MP2 correlation 
energy Based on the HF noninteracting Green function [equation (124]) with HF orbitals] , 
the 2nd-order self-energy can be expressed as 

£( 2 )(1,2)= -i J c/3c/4G HF (l,2)G HF (3,4)G HF (4,3)t;(l,3)t;(2,4) + 

% y"rf3rf4G HF (l,4)G HF (4,3)G HF (3,2)t;(l,3)t;(2,4). (31) 

Here for notational simplicity we have used 1 = (ri, tx), and v(l, 2) = v(ri — r 2 )5(ti, to)- 
After integrating out the internal (spatial and time) coordinates in ( 13T]) . the final 
expression for T,^ (in frequency domain) within the HF molecular orbital basis is given 
by 

occ unocc 
I a p,a' 



occ unocc 

-53 53 5^^ mp ' a i Za ' a )^' cr i an ' a 



IX 
a p 

9(e F - ep) 9{e p - e F ) 



(32) 



u + e au - ei a - e pa - vq u + e ta - e aa - e pr7 + ir]_ 

where 9{x) is the Heaviside step function and 77 — )- + . Again, the two lines in 
correspond to correlations arising from the 2nd-order direct and 2nd-order exchange 
interaction, respectively. With this self-energy single-particle excitation energies - here 
denoted MP2 quasiparticle energies - can be obtained by adding a correction to the HF 
orbital energies 



r • - - pHF 1 V(2) I MP2N /oo\ 



,MP2 
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Similar to the MP2 correlation energy, the quality of the MP2 self-energy as given by 
fl32|) will deteriorate as the single-particle energy gap shrinks. 

By means of the Dyson equation G = G° + G°SG self-consistent Green functions 
and self-energies could be obtained. The advantage is that the self-consistent Green 
function will then be independent of G° and satisfies particle number, momentum and 
energy conservation |122[ 1123] . For an in-depth discussion, we refer to |124[ 1125] and 
references therein. 

3. Resolution of identity for HF, MP2, RPA, and GW 

3.1. Background 

In this section we present the basic resolution of identity (RI) formalism: the auxiliary 
basis sets, different variants of RI, and the working equations for HF, hybrid density 
functionals, MP2, RPA, and GW. Similar accounts have been given in the literature 
on one or another of the above methods, and we encourage readers to consult them 
[761 EO CHI ESI 1271 [1281 129]. Our aim here is to lay out a complete description of 
all necessary specifics in one consistent notation that will allow us to present our own 
developments (next sections) in a self-contained way. 

The common ingredient of all techniques introduced in section [2] (Hartree-Fock, 
hybrid functionals, MP2, RPA, GW) are the four-orbital Coulomb integrals 

{mn, a\ab, a') = ^(ij l^c^c^c^c^,, (34) 

ijkl 

where (ij\kl) are the two-electron integrals in a basis set representation as defined in (I14p . 
and c l ma are the eigen-coefficients for the molecular orbitals. Computing (and possibly 
storing) these 4-center 2-electron integrals can be a major bottleneck for approaches 
beyond LDA and GGAs. 

For analytical GTOs algorithms have been developed to handle (mn, a\ab, a') 
integrals efficiently and on-the-fly |130[ 1131] . More general NAOs, however, are not 
amenable to such algorithms. In the context of HF, we note that RI is not the only 
technique available to deal with the four-center integrals: Making use of the translational 
properties of spherical harmonics, Talman and others |132[ I133[ I134[ 1135] have developed 
techniques based on multipole expansions of basis functions. Multi-center NAO integrals 
can then be treated partially analytically. Alternatively, efficient Poisson solvers [136J 
have recently been used to enable direct NAO HF calculations through four-center 
integrals for simple systems |133[ 1136] . Finally, a yet different numerical route (based on 
expanding orbital products directly) has been adopted along the line of time- dependent 
DFT in the linear- response framework |137] . RI is, however, most successful at reducing 
the computational load compared to direct four-center integral based methods, most 
prominently for MP2 in quantum chemistry [77], which is why we pursue this route 
here. 
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3.2. Auxiliary basis 

The resolution of identity (RI) or synonymously density fitting technique [Til E21 EH1 EU 
EE [76j [TTJ EH] amounts to representing pair products of atomic basis functions <fi(r)(pj(r) 
in terms of auxiliary basis functions (ABFs), 

Pij (r) = ^(r)^(r) « p^(r) = ^TC^(r). (35) 

/x = 1, 2, . . . , iVaux labels the auxiliary basis functions {P^}, are the expansion 
coefhcents, and P ij(r) and pij{v) here denote pair products of basis functions and their 
approximate expansion in ABFs. The evaluation of the 4-center integrals in ([H]) then 
reduces to 

(y|W)» J2 C ?M») C m> ( 36 ) 

(A*|^) =V IW = f Zp^P- drdr'. (37) 
y | r — r | 

To determine the expansion coefficients C^, three-center integrals involving the ABFs 
and the pair products of the NAOs are required. Thus the expensive (both in time and 
memory, if there is a need to pre-compute numerical matrix elements) 4-center integrals 
reduce to the much cheaper 3-center and 2-center ones in RI. The key reason for the 
success of RI lies in the fact that the set of all possible pair products {y?i(r) ■ ^ (r)}, as a 
set of basis functions in three-dimensional function space, is heavily linearly dependent. 
Their number scales quadratically with system size, while a non-redundant basis set 
that expands the same three-dimensional space should scale linearly with system size. 
For example, the non-interacting response function x° i n f[19l) - as well as the screened 
Coulomb interaction W in f [26l) . is written in terms of orbital pair products, and hence 
can be represented in terms of the ABFs. As will be shown below, this naturally leads 
to a RI implementation for RPA and GW . 

Next we will present RI formulations for all pertinent methods in this paper before 
presenting our specific choice for the ABFs subsequently. 

3. 3. Metric and variational principle in RI 

For a given set of ABFs {P M (r)}, the way to determine the expansion coefficients 
is not unique. Different variational procedures give rise to different versions of RI and 
different working equations for computing the [7U [721 E3 EH E2, [761 E3 EH] • 
The expansion error of basis products in terms of the ABFs [equation f[35l) ] is 

5 Pij (r) = ~ Pij (r) - Pij (r) = £ C^(r) - ^(r)^(r) . (38) 

One choice for the construction of the expansion coefficients is to minimize the 
residual 5 P ij(r). A simple least-square fit amounts to minimizing the norm of the 
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residual, J \6pij(r)\ 2 dr and yields 

V 

where 

(ij» = / Vi {v) Vj {v)PJv)dv , (40) 



and 

S M „ = / P^{v)P v {v)dv. (41) 



Combining (j36J) and fl39|) . one arrives at the following approximation to the four-center 
Coulomb integrals, 

(ij\kl) = J2 (i3\»)S^V w 'S;,l,W) . (42) 

In the literature equation (142 j) is therefore referred to as the "SVS" version [75] of RI 
( "RI-SVS" in the following) because of the appearance of the inverse S before and after 
the V matrix in (j4"2"j) . 

A better criterion for obtaining C£ is to directly minimize the RI error of the 
4-center integrals themselves, 

$Iij,ki = (PijlPki) - (Pij\Pki) ■ (43) 
As shown by Whitten [72] has an upper bound: 



< (Spij 1 5 pa ) 1/2 (8p kl 1 8p kl ) 1/2 + (p w 1 5p k i) 1/2 (fy f , | <Jp y ) 1/2 + ( py | Spij ) 1/2 (<Jp w 1 5p kl ) 1/2 . 

(44) 

The minimization of W^jy can thus be achieved by independently minimizing 5Uij = 
(5pij\6pij) and SUki = (Spki\$Pki)- the self- repulsion of the basis pair density residuals. 
Minimizing SUy with respect to C£ leads to [721 1731 1741175] 

q} = £(vW?, ( 45 ) 

V 

where 

(m = jj mMwyi dIdI ~ , (46) 

and the Coulomb matrix V is defined in (13"71) . Combining (1361) and (1451) one obtains the 
following decomposition of the 4-center ERIs 

(y|*0 » ■ (47) 

Equation (T47J is based on the global Coulomb metric and corresponds to the "RI-V" 
method [75j mentioned earlier. "RI-V" has long been known to be superior to "RI-SVS" 
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in the quantum chemistry community [721 [TjH EH ES] • This can be most easily seen by 
inspecting the error introduced in the self- Coulomb-repulsion of the NAO pairs, 

= (Pij\Pij) ~ (PijlPij) = 2($Pij\Pij) ~ (SPij\$Pij) ■ ( 48 ) 
In the RI-V approximation, the first term in the above equation vanishes, and the non- 
zero contribution comes only from the second order of 5pij. This can be readily verified 
as follows, 

(&Pij\Pij) = (PijlPij) ~ (PijlPij) = J2 C ii V ^ C ij ~ ^2^W) C ij 

= £(v»c& - £(v» c s = 

V V 

where ( )35l) and particularly (145 p have been used. In contrast, in RI-SVS the term linear 
in dp is non-zero and represents the dominating contribution to the total error. 

Our preferred flavor of RI is therefore RI-V, based on the Coulomb metric 
[721 1731 1741 1751 1761 1771 178] , on which all working equations for HF and other approaches 
presented further down in this section are based. Before proceeding we reiterate that 
RI-V continues to be the de facto standard in quantum chemical calculations, due to 
its well-established accuracy and reliability J77J [78] . 

That said, the long-range nature of the Coulomb interaction does present a 
bottleneck for implementations that scale better than the textbook standard (e.g., 
better than 0(N 4 ) for Hartree-Fock). In order to avoid delocalizing each localized 
two-center basis function product entirely across the system through C^, more localized 
approaches would be desirable. Research into better-scaling RI expansions that retain 
at least most of the accuracy of the Coulomb metric is thus an active field, for example 
by Cholesky decomposition techniques or an explicitly local treatment of the expansion 
of each product (see, e.g., [1381 EH HH ED EM EM HO US] for details). In fact, a 
promising Coulomb-metric based, yet localized, variant of RI has been implemented in 
FHI-aims. In this approach products of orbital basis functions are only expanded into 
auxiliary basis functions centered at the two atoms at which the orbital basis functions 
are centered, but the appropriate RI sub-matrices are still treated by the Coulomb 
metric [146J. As expected, the error cancellation in this approach is not as good as 
that in full RI-V, but — for Hartree-Fock and hybrid functionals — certainly more than 
an order of magnitude better than in RI-SVS, creating a competitive alternative for 
cases where RI-V is prohibitive. More details would go beyond the scope of this paper 
and will be presented in a forthcoming publication [147] . 

3.4- HF and hybrid functionals 

The key quantity for HF and hybrid functionals is the exact-exchange matrix - the 
representation of the non-local exact-exchange potential [equation ([8])] in terms of basis 
functions as given in [T2l). Its RI-V expansion follows by inserting ( 1471) into ( TT2l : 

kl puV fi kl 
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where 

= Ew^ 2 = E c ^ ■ ( 5 °) 

V V 

£* )tT must thus be recomputed for each iteration within the self-consistent field (scf) 
loop. The required floating point operations scale as ■ iV aux . The transformation 
matrix M^ k , on the other hand, is constructed only once (prior to the scf loop), requiring 
N£ ■ N^ ux operations. 

The numerical efficiency can be further improved by inserting the expression for 
the density matrix (ITS]) into (149]) : 



occ: 



= E E E E = E E B ^ B U ■ (si) 

n n \ k / \ I / n 11 

The formal scaling in (151 p is now iV occ • N-^ 2 ■ N aux , with iV occ being the number of 
occupied orbitals, and thus improved by a factor N h /N occ (typically 5 to 10). Once the 
exact-exchange matrix is obtained, the exact-exchange energy follows through 

E ^ = -\Y,^ D ^- ( 52 ) 



For a variety of physical problems, combinding HF exchange with semi-local 
exchange and correlation of the GGA type gives much better results than with pure 
HF or pure GGAs [29J. Various flavors of these so-called hybrid functionals exist in the 
literature. The simplest one-parameter functionals are of the following form 

ET = £ X G C GA + «(^ HF - £ X GGA ) . (53) 

In the PBEO hybrid functional [30], the GGA is taken to be PBE, and the mixing 
parameter a is set to 1/4. Naturally, the computational cost of hybrid functionals is 
dominated by the HF exchange. Once HF exchange is implemented, it is straightforward 
to also perform hybrid functional calculations. 



3.5. MP2 (total-energy correction and self-energy) 

To compute the MP2 correlation energy in (FToT) and the MP2 self-energy in (132]) using 
the RI technique, the MO-based 4-orbital 2-electron Coulomb integrals are decomposed 
as follows 

(ma,a\nb,a') = J2°^,A !al . (54) 
The 3-orbital integrals can be evaluated by 

OL,. = E M £ c -cL- (55) 

ij 



18 



following (El]), ( HTj) . and (l5T?j) . Plugging (1511) into (ITS]) , one obtains the RI-V version of 
the MP2 correlation energy 

1 occ unocc / \ 

mn ab 0,0' V A 4 / 

(56) 

In practice, one first transforms the atomic orbital-based integrals to MO-based 
ones 0^ naa . The transformation scales formally as iV occ • N£ ■ iV aux and the summation 
in fl56|) as N 2 CC • (N h — N occ ) 2 ■ iV aux . Like in Hartree-Fock, the scaling exponent 0(N 5 ) is 
therefore not reduced by RI-V. However, the prefactor in RI-MP2 is one to two orders 
of magnitude smaller than in full MP2. 

The computation of the MP2 self-energy at each frequency point proceeds 

analogously to that of the correlation energy. In our implementaton, we first calculate 

(2) 

the MP2 self-energy on the imaginary frequency axis, S^ io ., and then continue 
analytically to the real frequency axis using either a "two-pole" model |119j . or the 
Pade approximation. Both approaches have been implemented in FHI-aims and can 
used to cross-check each other to guarantee the reliability of the final results. 

3.6. RPA and GW 

To derive the working equations for the RPA correlation energy and the GW self-energy 
in the RI approximation, it is illuminating to first consider the Rl-decomposition of 
Tr[x° \iu)v\. Combining ( !T9|) and (|54|) we obtain 

occ unoc c 

Tr [ x vh = E E E E Yu -T Te + cx - ' (57) 

^ tacr T t mcr 

fi a m a 

where 0^ naa is given by (1551) . Next we introduce an auxiliary quantity U(iu): 



occ unocc 



n(- v = E E E °t°°l m e + c - c - ' (K 

tacr T t mcr 



a m a 



which allows us to write 

Tr [x°(iu)v] = Tr [v l ' 2 X \iu)v 1/2 ] = Tr [U(iu)] . (59) 

Thus the matrix II can be regarded as the matrix representation of the composite 
quantity v 1 ' 2 x°{iw)v 1 ' 2 using the ABFs. It is then easy to see that the RPA correlation 
energy ( ITHj) can be computed as 



1 r°° 

e rpa = ± / duTr p n (1 _ u ^ + n ^ 

27T Jo 
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1 f°° 

= — duj {In [det(l - n(iw)) + Tr [U(iu)}} . 
^ Jo 

(60) 

using the general property Tr[ln(A)] = In [det(A)] for any matrix A. This is very 
convenient since all matrix operations in |60] occur within the compressed space of ABFs 
and the computational effort is therefore significantly reduced. 

In practice, we first construct the auxiliary quantity U(ico) [equation ( 1581) ] on a 
suitable imaginary frequency grid where we use a modified Gauss- Legendre grid (see 
appendix [7] for further details) with typically 20-40 frequency points. For fixed frequency 
grid size, the number of required operations is proportional to N occ • N unocc ■ N^ ux (iV unocc 
is the number of unoccupied orbitals using the full spectrum of our Hamiltonian matrix). 
The next step is to compute the determinant of the matrix 1 — H(iu) as well as the trace 
of U(iu>). What remains is a simple integration over the imaginary frequency axis. Thus 
our RI-RPA implementation is dominated by the step in ( )58l) that has a formal 0(N 4 ) 
scaling. An (9(iV 4 )-scaling algorithm of RI-RPA was recently derived from a different 
perspective |129j , based on the plasmonic formulation of RPA correlation energy |107] 
and a transformation analogous to the Casimir- Polder integral |148j . An even better 
scaling can be achieved by taking advantage of the sparsity of the matrices involved 

Finally we come to the RI-V formalism for GW. To make ( 130]) tractable, we 
expand the screened Coulomb interaction W(iu) in terms of the ABFs. Using (12?]) and 
U(ico) = v 1 / 2 Xo(^)'^ 1//2 , we obtain 

W^(iu) = Jj drdr'P;(r)W(r, r', iw)P„(r) = E V$[l - U(iu)]^vX 2 ■ 

To apply the Rl-decomposition to ( 1301) . we expand 

C«Vw(r) = E E ^( f )^C4 , (61) 

ij A 4 

where (ITOjl and (J35l) are used. Combing (!30l) . ( )50l) . ( |55l) . and El]) gives 

(nm, o\W\iu)\mn, a) = E E °w t 1 " ■ (62) 

Inserting ( 1621) into (|29l) . one finally arrives at the RI- version of the G°W° self-energy 

"OO j 

" li.r)= N / dr.i' 

na 



-y r 



^: w °uuj) = - - > - / du>, — — x 



27r J^oo ibJ + iuj' + cf — e r 

E^^nNCoL,,. (63) 

As stated above for the MP2 self-energy, the expression is analytically continued to the 
real- frequency axis before the quasiparticle energies are computed by means of ([28]) . 
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4. Atom-centered auxiliary basis for all-electron NAO calculations 



4-1. Orbital basis set definitions 

For the practical implementation, all the aforementioned objects (wave functions, 
effective single-particle orbitals, Green function, response function, screened Coulomb 
interaction etc.) are expanded either in a single-particle basis set or an auxiliary basis 
set. We first summarize the nomenclature used for our NAO basis sets [64] before 
defining a suitable auxiliary basis prescription for RI in the next subsections. 

NAO basis sets {^(r)} to expand the single-particle spin orbitals ip n a{ r ) [equation 
( jTUl) ] are of the general form 



u s { a )iK is a radial function centered at atom a, and YJ m (f a ) is a spherical harmonic. The 
index s(a) denotes the element species s for an atom a, and k enumerates the different 
radial functions for a given species s and an angular momentum I. The unit vector 
r a = (r — R a )/|r — R a | refers to the position R a of atom a. The basis index i thus 
combines a, k, I, and m. 

For numerical convenience, and without losing generality, we use real- valued basis 
functions, meaning that the Yi m (Q) denote the real (for m — 0, • • • , I) and the imaginary 
part (for m = —I, - ■ ■ , — 1) of complex spherical harmonics. For NAOs, u s ^i K (r) need 
not adhere to any particular analytic shape, but are tabulated functions (in practice, 
tabulated on a dense logarithmic grid and evaluated in between by cubic splines). Of 
course, Gaussian, Slater-type, or even muffin-tin orbital basis sets are all special cases 
of the generic shape (l64"j) . All algorithms in this paper could be used for them. In 
fact, we employ the Dunning GTO basis sets (see [149\ 1150] and references therein) for 
comparison throughout this work. 

Our own implementation, FHI-aims, [64] provides hierarchical sets of all-electron 
NAO basis functions. The hierarchy starts from the minimal basis composed of the 
radial functions for all core and valence electrons of the free atoms. Additional groups 
of basis functions, which we call tiers (quality levels) can be added for increasing 
accuracy (for brevity, the notation is minimal, tier 1, tier 2, etc.). Each higher level 
includes the lower level. In practice, this hierarchy defines a recipe for systematic, 
variational convergence down to meV/atom accuracy for total energies in LDA and 
GGA calculations. The minimal basis (atomic core and valence radial functions) is 
different for different functionals, but one could as well use, e.g., LDA minimal basis 
functions for calculations using other functionals in cases their minimal basis functions 
are not available. We discuss this possibility for HF below. 

To give one specific example, consider the nitrogen atom (this case and more are 
spelled out in detail in Table 1 of Ref. [64J). There are 5 minimal basis functions, 
of Is, 2s, and 2p orbital character, respectively. In a shorthand notation, we denote 
the number of radial functions for given angular momenta as (2s lp) (two s-type radial 




(64) 



r 
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functions, one p-type radial function). At the tier 1 basis level, one s, p, and d function 
is added to give a total of 14 basis functions (3s2pld). There are 39 basis functions 
(4s3p2dlflg) in tier 2, 55 (5sAp3d2flg) in tier 3, and 80 (6s5p4d3f2g) in tier 4. 

4-2. Construction of the auxiliary basis 

In the past, different communities have adopted different strategies for building auxiliary 
basis sets. For the GTO-based RI-MP2 method [77], which is widely used in the 
quantum chemistry community, a variational procedure has been used to generate 
optimal gaussian-type atom-centered ABF sets. In the condensed matter community 
a so-called "product basis" has been employed in the context of all-electron GW 
implementations based on the linearized muffin-tin orbital (LMTO) and/or augmented 
plane- wave (LAPW) method [15 lj to represent the response function and the Coulomb 
potential within the muffin-tin spheres |126j 11521 1127j . Finally, it is even possible to 
generate ABFs only implicitly, by identifying the "dominant directions" in the orbital 
product space through singular value decomposition (SVD) [144[ I145[ 1153] . As will be 
illustrated below, our procedure to construct the ABFs combines features from both 
communities. Formally it is similarly to the "product basis" construction in the GW 
community, but instead of the simple overlap metric the Coulomb metric is used to 
remove the linear dependence of the "products" of the single-particle basis functions, 
and to build all the matrices required for the electronic structure schemes in this paper. 

Our procedure employs numeric atom-centered ABFs whereby the infrastructure 
that is already available to treat the NAO orbital basis sets can be utilized in many 
respects. Specifically the ABFs are chosen as 

P„(r) = ^^Y lm (h) (65) 

just like for the one-particle NAO basis functions in ( 1641) . but of course with different 
radial functions. To distinguish the auxiliary basis functions from the NAO basis 
functions we denote the radial functions of the ABFs as £ s (a)i K - 

The auxiliary basis should primarily expand products of basis functions centered on 
the same atom exactly, but at the same time be sufficiently flexible to expand all other 
two-center basis function products with a negligible error. In contrast to the ABFs used 
in the GTO-based RI-MP2 method [77] , in our case the construction of auxiliary basis 
functions follows from the definition of the orbital basis set. At each level of NAO basis, 
one can generate a corresponding ABF set, denoted as aux-min, aux-tier 1, auxdier 2, 
etc. We achieve this objective as follows: 

(i) For each atomic species (element) s, and for each I below a limit Z™ ax , we 
form all possible "on-site" pair products of atomic radial functions {£ s z K ( r ) = 
Uskxh ( r ) u sk 2 h( r )} ■ The allowed values of I are given by the possible multiples of 
the spherical harmonics associated with the orbital basis functions corresponding 
to u sklh and u sk2h) i.e., \l x - l 2 \ < I < \h + h\- 
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(ii) Even for relatively small orbital basis sets, the number of resulting auxiliary radial 
functions is large. They are non-orthogonal and heavily linear dependent. 
We can thus use a Gram-Schmidt like procedure (separately for each s and I) to 
keep only radial function components £ S / K (r) that are not essentially represented 
by others, with a threshold for the remaining norm e orth , below which a given 
radial function can be filtered out. In doing so the Coulomb metric is used in 
the orthogonalization procedure. The result is a much smaller set of linearly 
independent, orthonormalized radial functions {£sz K (V)} that expand the required 
function space. 

(iii) The radial functions {£,s(ci),Ik} are multiplied with the spherical harmonics Yi m (r a ) 
as in fl65l) . 

(iv) The resulting {P M (r)} are orthonormal if they are centered on the same atom, but 
not if situated on different atoms. Since we use large ABF sets, linear dependencies 
could also arise between different atomic centers, allowing us to further reduce the 
ABF space through SVD of the applied metric (S in the case of RI-SVS, V in the 
case of RI-V). For the molecule- wide SVD we use a second threshold e svd , which is 
not the same as the on-site Gram-Schmidt threshold e orth . 

For a given set of NAOs, the number of the corresponding ABFs depends on the angular 
momentum limit /™ ax in step 1 and the Gram-Schmidt orthonormalization threshold 
£ orth , and to a small extent on e svd . For RI-V, as documented in the literature [78J and 
demonstrated later in this work (section l4.5p . it is sufficient to keep Z™ ax just one higher 
than the highest angular momentum of the one-electron NAOs. Usually £ orth = 10~ 2 or 
10 -3 suffices for calculations of energy differences. Nevertheless both /™ ax and e orth can 
be treated as explicit convergence parameters if needed. In practice, we keep e svd as 
small as possible, typically 10 -4 or 10~ 5 , only large enough to guarantee the absence of 
numerical instabilities through an ill-conditioned auxiliary basis. The resulting auxiliary 
basis size is typically 3-6 times that of the NAO basis. This is still a considerable size and 
could be reduced by introducing optimized ABF sets as is sometimes done for GTOs. 
On the other hand, it is the size and quality of our auxiliary basis that guarantees 
low expansion errors for RI-V, as we will show in our benchmark calculations below. 
We therefore prefer to keep the safety margins of our ABFs to minimize the expansion 
errors, bearing in mind that the regular orbital basis introduces expansion errors that 
are always present. 

4-3. Numerical integral evaluation 

With a prescription to construct ABFs at hand, we need to compute the overlap integrals 
C^, defined in fl39|) for RI-SVS or in f )45l) for RI-V, respectively. We also need their 
Coulomb matrix given by fl37|) in general, and additionally the "normal" overlap matrix 
S^ v given by (jHJ) in RI-SVS. Having efficient algorithms for these tasks is enormously 
important, but since many pieces of our eventual implementation exist in the literature, 
we here only give a brief summary and refer to separate appendices for details. 
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Since our auxiliary basis set {P^} is atom-centered, we obtain the Coulomb 
potential Q M (r) of each P^(r) by a one-dimensional integration for a single multipole 
(J7j). The required three-center integrals 



are carried out by standard overlapping atom-centered grids as used in many quantum- 
chemical codes for the exchange-correlation matrix in DFT [154[ l65 | l64" tll55| , seed The 
same strategy works for two-center integrals 



As an alternative, we have also implemented two-center integrals following the ideas 
developed by Talman |133[ 1134] , which are described in [7J and UJ In summary, we thus 
have accurate matrix elements at hand that are used for the remainder of this work. 

4-4- Accuracy of the auxiliary basis: expansion of a single product 

In this section we examine the quality of our prescription for generating the ABFs as 
described in section 14.21 Our procedure guarantees that the ABFs accurately represent 
the "on-site" products of the NAO basis pairs by construction, but it is not a priori clear 
how the "off-site" pairs are represented. The purpose of this section is to demonstrate 
the quality of our ABFs to represent the "off-site" pairs. 

In the left panels of figure [H we plot p2s-2 Px (r) for a simple N 2 molecule at the 
equilibrium bonding distance (d = 1.1 A) - the product of the atomic 2s function 
from the left atom and the atomic 2p x function from the right atom. We compare this 
directly taken product to its ABF expansions, both in RI-SVS [equation (|39|) ] and in 
RI-V [equation (|45l) ]. The particular product P2s-2p x { r ) ls P ar ^ °f the minimal basis 
of free-atom like valence radial functions. As we increase the orbital basis set (adding 
tier 1, tier 2, etc.), the exact product remains the same, whereas its ABF expansion 
will successively improve, since the auxiliary basis set is implicitly defined through 
the underlying orbital basis. Three different levels of ABF sets are shown (aux-min, 
aux-tier 1, auxMer 2 from the top to bottom panels). The onsite threshold e orth is 
set to 10" 2 , and the global SVD threshold e svd is set to 10" 4 , yielding 28, 133, and 
355 ABFs, respectively. In the right panels of figure [fl the corresponding Sp 2 s-2 Px ( r ) ~ 
the deviations of the ABF expansions from the reference curve - are plotted. One can 
clearly see two trends: First, the quality of the ABF expansion improves as the number 
of ABFs increases. Second, at the same level of ABF, the absolute deviation of the RI-V 
expansion is larger than the RI-SVS expansion. This is an expected behaviour for the 
simple pair product: RI-V is designed to minimize the error of the Coulomb integral, 
see section 13.31 In either method, the remaining expansion errors are centered around 
the nucleus, leading to a relatively small error in overall energies (in the 3-dimensional 
integrations, the integration weight r 2 dr is small). 
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Figure 1. (Color online) Left panels: the product of the atomic 2s and 2p x 
functions centered respectively on the two atoms in a N2 molecule along the bonding 
direction, and its approximate behaviours from the RI-SVS and RI-V expansions for 
three hierarchical levels of ABFs (from auxjmin to aux-tier 2). Right panels: The 
corresponding deviations of the RI-SVS and RI-V expansions from the reference curve. 
The positions of the atoms are marked by blue dots at the cc-axis. 

Table [1] gives the errors of the Coulomb repulsion Slij^j of the 2s- 2p x NAO basis 
pair under the RI approximation for the three ABF basis sets of figure [TJ Table [TJ also 
includes the influence of the threshold parameters e orth and e svd , separate for RI-SVS 
(top half) and RI-V (bottom half). The error diminishes quickly with increasing ABF 
basis size, and is 2-3 orders of magnitude smaller in RI-V than in RI-SVS. By decreasing 
£ ort h, the number of ABFs at each level increases, improving particularly the accuracies 
at the auxjmin and auxJier 1 levels. The global SVD threshold e svd comes into play 
only for the larger basis sets. In general, both control parameters have a much bigger 
effect on RI-SVS than RI-V, underscoring the desired variational properties of RI-V 
[H1E21II31I711E5]. 

4-5. Accuracy of the auxiliary basis: energies and thresholds 

Next we turn to the accuracy of our auxiliary basis prescription for actual HF and MP2 
(total and binding energy) calculations. For this purpose, we employ all-electron GTO 
basis sets, when possible, to be able to refer to completely independent and accurate 
implementations from quantum chemistry without invoking the RI approximation 
(referred to as "Rl-free" in the following). Our specific choice here is the GTO-based 
NWChem [156J code package, where "Rl-free" results can be obtained using traditional 
methods of quantum chemistry. In the following we compare our Rl-based HF and MP2 
results with their "Rl-free" counterparts produced by NWChem, in order to benchmark 
the accuracy of the RI implementation in FHI-aims. All results presented in this section 
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Table 1. The errors S l2s-2p m ,2s-2p^ 

(in eV) introduced in the RI-SVS and RI-V for 
calculating the self-repulsion of NAO pair products (p2s-2p :c |/°2s-2p x 

) for N 2 (d = 1.1 A) 

at three levels of ABF basis sets. The number of ABFs that survive the SVD is also 
shown. 
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correspond to the cc-pVQZ basis set of Dunning et al. |157[ I149[ 1150] The convergence 
behaviour with respect to NAO / GTO single-particle basis set size will be the topic of 
next section. 

We first check the quality of the ABF prescription for light (first and second- 
row) elements. We again choose N2 as a first illustrating example. Table [2] presents 
RI-HF total and binding energy errors for the N2 molecule at bond length d=l.l A. 
The reference numbers given at the bottom of the table are from "Rl-free" NWChem 
calculations. All other numbers were obtained with FHI-aims and the ABF prescription 
of section 14.21 Total and binding energy errors are given for several different choices 
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Table 2. The deviation of HF total (A£" t ot) and binding energies (AE D ) (in meV) 
from NWChem reference values for N 2 at bond length d =1.1 A. RI-V and RI-SVS 
calculations are done using the Dunning cc-pVQZ basis set. |157l I149j The reference 
E tot and E^ values (in eV) from NWChem calculations are shown at the bottom. The 
binding energies are BSSE-corrected using the counterpoise method. 
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of thresholding parameters e orth and e svd (see section H~2l) . All binding energy errors 
are obtained after a counterpoise correction [ST] to remove any possible basis set 
superposition errors (BSSE) (see also lowest panel of figure [2] below) . 

For N2 at equilibrium bonding distance, the salient results can be summarized as 
follows. First, we see that RI-V with our ABF prescription implies total energy errors 
for Hartree-Fock of only ~0. 1-0.2 meV, while the corresponding RI-SVS errors are much 
larger. Both methods can, however, be adjusted to yield sub- me V binding energy errors, 
with those from RI-V being essentially zero. This is consistent with our observations for 
the error in the self-repulsion integrals in section 14.41 Since RI-V performs much better 
than RI-SVS, we will only report RI-V results for the remainder of this paper. 

The excellent quality of our ABFs is not restricted to the equilibrium region of N 2 . 
In the top panel of figure [2] the restricted HF total energies are plotted for a range of 
bonding distances. The RI-V numbers are in very good agreement with the reference 
throughout. For greater clarity, the total-energy deviation of the RI-V result from 
the reference is plotted in the middle panel of figure [2j One can see that the total- 
energy errors are in general quite small, but the actual sign and magnitude of the errors 
vary as a function of bond length. The deviation is shown for two different choices of 
the ABFs, (e orth =10~ 2 , e svd =10~ 4 ) (standard accuracy) and (e OTth = KT 3 , £ svd =lCr 5 ) 
(somewhat tighter accuracy). It is evident that the tighter settings produce a smoother 
total energy error at the sub-meV level, but, strikingly, there is no meaningful difference 
for the counterpoise corrected binding energy of N2 (bottom panel of figure [2]). 

Next we demonstrate the quality of our ABFs for a set of molecules consisting 
of first and second-row elements. In figure [3] and H] the non-relativistic RI-V HF and 
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Figure 2. (color online) Upper panel: RI-V HF total energies as a function of bond 
length for N 2 , in comparison with NWChem reference values. Middle panel: the 
deviation of the RI-V HF total energies from the reference values for two sets of 
thresholding parameters. Lower panel: the deviation of the BSSE-corrected RFV HF 
binding energies from the reference values for the same sets of thresholding parameters. 
The cc-pVQZ basis is used in all the calculations. Note that the dependence of the 
total energies on the thresholding parameters is not visible in the upper panel. 




Figure 3. (color online) Deviations of RFV HF total energies (red circles) and 
atomization energies (blue squares) from the corresponding reference values for 20 
small molecules. The three panels illustrate the dependence of the RI errors on the 
truncation parameters e orth , e svd , and the highest ABF angular momentum l ABF ~ max 
^AO-max denotes the highest angular momentum of the single-particle atomic orbitals). 
Experimental equilibrium geometries and the gaussian cc-pVQZ basis are used. 
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Figure 4. (color online) Deviations of RI-V MP2 total energies and atomization 
energies from the corresponding reference values for 20 small molecules. Nomenclature 
and labelling are the same as figure |3] 



MP2 total energy errors and atomization energy errors for 20 molecules are shown. In 
all cases the total energy error is below 1 meV/atom, demonstrating that the meV- 
accuracy in total energy can routinely be achieved for the RI-V approximation with 
our ABFs for light elements. In addition, it is evident that varying the auxiliary 
basis convergence settings has essentially no influence on the low overall residual error, 
which is attributed to other small numerical differences between two completely different 
codes (analytical integrations in NWChem vs. numerical integrations in FHI-aims, for 
example). Furthermore, it is also clear that it is sufficient to choose the highest angular 
momentum in the ABF construction (/ ABF - max ) to be just 1 higher than that of the 
single-particle atomic orbitals (/ A °- max ). 

Having established the quality of our ABFs for light elements, we now proceed 
to check their performance for the heavier elements where some noteworthy feature is 
emerging. In figure [5] we plot the errors in the RI-V HF total energies and binding 
energies for Cu 2 as a function of the bond length. Again the cc-pVQZ basis and the 
NWChem reference values are used here. Using the thresholding parameters (£ orth =10~ 2 , 
£ svd =10 -4 ), the RI-V HF total energy error can be as large as ~ 15 meV/atom for 
copper, in contrast to the < 1 meV/ atom total energy accuracy for light elements. This 
is because for Cu, deep core electrons are present and the absolute total-energy scale is 1- 
2 orders of magnitude larger than that of light elements. The residual basis components 
eliminated in the on-site Gram-Schmidt orthonormalization procedure (see section 14.2ft 
thus give bigger contributions to the total energy on an absolute scale (although not 
on a relative scale). Indeed by decreasing e orth the total-energy error gets increasingly 
smaller, and 1-1.5 meV/atom total-energy accuracy can be achieved at e orth = 10~ 4 , 
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Figure 5. (color online) Upper panel: Deviations of non-relativistic RI-V HF total 
energies from the NWChem reference values as a function of bond length for Cu2 for 
four sets of thresholding parameters. Lower panel: the deviation of the BSSE-corrected 
RFV HF binding energies from the reference values for the same sets of thresholding 
parameters. The cc-pVQZ basis is used in all the calculations. 

as demonstrated in the upper panel of figure [H In contrast, similar to the case of 
light elements, the errors in the BSSE-corrected binding energies are significantly below 
0.1 meV/atom along a large range of bonding distances, regardless of the choice of 
thresholding parameters. And also increasing the highest angular momentum for ABFs 
beyond / A °- max 4- 1 does not give noticeable improvements. 

Finally we look at the quality of our ABFs for even heavier elements - the Au 
dimer. Due to the strongly localized core states in Au, all-electron GTO basis sets 
that are converged to the same level of accuracy as for N and Cu above are, to our 
knowledge, not available for Au. Furthermore, relativity is no longer negligible and 
must at least be treated at a scalar relativistic level. However, different flavors of 
relativistic implementations can differ heavily in their absolute total energy (even if 
all chemically relevant energy differences are the same). An independent reference for 
all-electron GTO-based HF total energy with the same relativistic treatments available 
in FHI-aims is not readily obtainable. Under such circumstances, we demonstrate here 
the total energy convergence with respect to our own set of thresholding parameters. 

In figure |6] we plot the deviations of the RI-V HF total and binding energies for AU2 
obtained with FHI-aims using NAO tier 2 basis with somewhat less tight thresholding 
parameters from those obtained with a very tight threshold setting (e orth = 10 -5 , 
£ svd _ T/he relativistic effect is treated using the scaled zeroth-order regular 

approximation (ZORA) |158j (see section H~5l) . but this detail is not really important for 
the discussion here. From figure [6] one can see that with thresholding parameters that 
are sufficient for light elements (e orth =10 _2 , £ svd =10 -4 ), the error in the total energy 
is even bigger - one order of magnitude larger than in the case of CU2. However, by 



30 



300 
200 
! 100 




o--o 
o--D--e i 



.-Q-O 

- B -B- -□- -Or 



o_ o_ . -o -o- o- « -o -o- <y $ 

A— A— A— A— A— A— A— A— A— A— A 



-« ❖ 



> 

s 





1 


_ — orth , „-2 svd , „-4 

o o e = io , e =10 


, 








0.2 




„ orth . _-2 svd , -5 

e-^e =10 , e =10 




1 1 


1 


0.1 




. . orth . _-3 svd , -5 

o-<>e =10 ,e =10 






a a orth , _-4 svd , „-5 

=10 ,e =10 




0.0 


§~S ~# r ff ~ - ■ a — S~ 










-0.1 




i,i 


, 




2 2.4 


2.8 3.2 
Bond length (A) 


3.6 4 



Figure 6. (color online) Upper panel: Deviations of scaled ZORA RI-V HF total 
energies from the reference values as a function of bond length for Au2 for four sets of 
thresholding parameters. The reference values here are also obtained with the RI-V 
HF approach with very tight thresholds (e orth = 10~ 5 ,£ svd = 10~ 5 ). Lower panel: the 
deviation of the BSSE-corrected RI-V HF binding energies from the reference values. 
The FHI-aims tier 2 basis and / ABF " max = jAO-max + 1 are uged in all the ca i cu i at i ns. 

going to tighter and tighter on-site ABF thresholding parameter e orth , the total-energy 
error can again to be reduced to the meV/atom level. Similar to the Cu2 case, the 
accuracy in the BSSE-corrected binding energy is still extraordinarily good, well below 
0.1 meV regardless of the choice of thresholding parameters. The counterpoise correction 
can thus be used, in general, as a simple, readily available convergence accelerator for 
binding energies. 

We conclude this section with the following remarks: our procedure for constructing 
the ABFs gives highly accurate and reliably results for the RI-V approximation. For light 
elements, one can readily get meV/ atom accuracy in total energies and sub-me V/ atom 
accuracy in binding energies, for a wide range of thresholding parameters. For heavy 
elements, meV/atom accuracy requires tigher thresholding parameters, particulary for 
the on-site orthonormalization e orth . However, sub-meV BSSE-corrected binding energy 
accuracy can always be obtained independent of the choice of thresholding parameters. 
Choosing converged yet efficient thresholding parameters thus obviously depends on the 
element in question. In fact, e OTth can be chosen differently for each element in the 
same calculation. For RI-V and light elements [Z = 1-10), we employ e ovth = 10~ 2 in 
the following. For heavier elements (Z > 18), we resort to e orth = 10 -4 , which yields 
negligibly small errors in total energies even for heavy elements, and e orth = 10 -3 for 
elements in between. The additional, system-wide SVD threshold e svd is set to 10~ 4 or 
tighter for the remainder of this paper. As shown above, its accuracy implications are 
then negligible as well. 
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5. NAO basis convergence for HF, hybrid density functionals, MP2, RPA, 
and GW 

Having established the quality of our ABFs for given orbital basis sets, we next 
turn to the quality of our actual NAO orbital basis sets for HF, hybrid density 
functionals, MP2, RPA, and GW calculations. Below we will separate the discussions 
of self-consistent ground-state calculations (HF and hybrid density functionals) and 
correlated calculations (MP2, RPA, and GW). As will be demonstrated below, with 
our basis prescription, a convergence of the absolute HF total energy to a high accuracy 
(mev/atom) is possible for light elements. In other cases (heavy elements or correlated 
methods), the total energy convergence is not achieved at the moment, but we show 
that energy differences, which are of more physical relevence, can be achieved to a high 
quality. 

5.1. HF and hybrid density functional calculations 

Here we will demonstrate how well the generic NAO basis set libraries described in 
section 14.11 and in [64J perform for HF and hybrid density functional calculations. As 
described earlier, our orbital basis sets contain a functional- dependent minimal basis, 
composed of the core and valence orbitals of the free atom, and additional functional- 
independent optimized basis sets {tiers). Thus, besides the discussion of the convergence 
behaviour of the generic optimized tiers basis sets, here we will also address the influence 
of the choice of the minimal basis which is in practice generated by certain atomic solver 
. For all-electron calculations for molecular systems with a given electronic-structure 
method, it would be best if the core basis functions were generated from the atomic solver 
using the same method. In this way, the behaviour of the molecular core wavefunctions 
in the vicinity of the nuclei would be accurately described at a low price. This is the 
case for LDA and most GGA calculations in FHI-aims. Similarly, for HF molecular 
calculations, it would be ideal if the minimal basis was generated by the HF atomic 
solver. Unfortunately at the moment the HF atomic solver is not yet available in our 
code, and in practice we resort to the minimal basis generated from other types of 
atomic solvers. This is not a fundamental problem, and the only price one has to pay 
is that more additional tiers basis functions are needed to achieve a given level of basis 
convergence. Nevertheless one should keep in mind that there is a better strategy here 
and in principle our basis prescription should work even better than what is reported 
here. 

In the following the NAO basis convergence for HF will be examined, and along the 
way the influence of the minimal basis will be illustrated by comparing those generated 
by DFT-LDA and Krieger-Li-Iafrate (KLI) [159j atomic solvers. The KLI method 
solves approximately the exact-exchange optimized-effective-potential (OEP) equation 
[160[ 1102] . by replacing the orbital- dependent denominator in the Green function (at 
zero frequency) appearing in the OEP equation by an orbital-independent parameter. 
This in practice reduces the computational efforts considerably without losing much 
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accuracy [161] . The KLI atomic core wavefunctions resemble the HF ones much better 
than the LDA ones do. As demonstrated below, by moving from the LDA minimal basis 
to the KLI ones in HF calculations, the abovementioned problem is alleviated to some 
extent. 

5.1.1. Light elements We first check the convergence behaviour of our NAO basis sets 
for light elements in HF and hybrid functional calculations. In figure [7] (left panel) the 
HF total energy of N 2 as a function of increasing basis set size is plotted, starting with 
two different sets of minimal bases - generated respectively from LDA and KLI atomic 
solvers. All other basis functions beyond the minimal part (the tiers) are the same 
for both curves. For comparison, the convergence behaviour of the LDA total energy 
of N2 is shown on the right panel for the same basis sets. The horizontal (dotted) 
lines indicate independently computed GTO reference values using NWChem and the 
Dunning cc-pV6Z basis set, which gives the best estimate for the HF total energy at 
the complete basis set limit |162[ 1163] . 

First, with both types of minimal basis the HF total energy can be systematically 
converged to within a few meV of the independent GTO reference value. This is 
reassuring, as we can thus use our standard NAO basis sets in a transferable manner even 
between functionals that are as different as LDA and HF. Furthermore, it is evident that 
the KLI-derived minimal basis performs better for HF, and similarly the LDA derived 
minimal basis performs better for LDA. As mentioned above, this is because the closer 
the starting atomic core basis functions to the final molecular core orbitals, the faster the 
overall basis convergence is. If the true HF minimal basis was used, we should expect 
an even faster basis convergence of the HF total energy, similar to the LDA total-energy 
convergence behaviour starting with the LDA minimal basis ((blue) circles in the right 
panel of figure Cj). In this case the BSSE in a diatomic molecular calculation should also 
be vanishingly small since the atomic reference is already accurately converged from the 
outset with the minimal basis. In practice, the reliance on KLI-derived minimal basis 
functions instead leads to some small BSSE-type errors in energy differences, as shown 
below. 

In figure M we present the NAO basis convergence of the HF binding energy for N2 
as a function of bond distance. Here we start with the KLI minimal basis and then 
systematically add basis functions from tier 1 to tier 4. Results are shown both without 
(left) and with (right) a counterpoise correction. A substantial improvement of the 
binding energy is seen between the tier 1 and the tier 2 basis, with only slight changes 
beyond tier 2. In the absence of BSSE corrections a slight overbinding is observed for 
tier 2 and tier 3. As noted above, we attribute this to the fact that in our calculations 
we used KLI core basis functions as a practical compromise and the atomic reference 
calculation is not sufficiently converged for the outset. The basis functions from the 
neighboring atom will then still contribute to the atomic total energy and this leads 
to non-zero BSSE. The counterpoise correction will cancel this contribution — which is 
very similar for the free atom and for the molecule — almost exactly. The counterpoise- 



33 



-2965.2 
-2965.3 
-2965.4 



u -2965.6 



o 
H 



-2965.7 
-2965.8 
-2965.9 







o-o KLI atomic solver 
O-G LDA atomic solver 




50 



100 




150 







50 



100 150 



Basis size 



Figure 7. (color online) NAO basis convergence test: HF and LDA total energies 
for N2 at d = 1.1 A as a function of increasing NAO basis set size (tier 1 to tier 4). 
The two convergence curves correspond to starting minimal basis sets generated by 
LDA and KLI atomic solvers, respectively. The dotted horizontal line marks the HF 
and LDA total energy computed using NWChem and the cc-pV6Z basis, which gives 
a reliable estimate of the basis-set limit within 2 mcV [162 . 
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Figure 8. (color online) HF binding energy of N2 as a function of the bond length for 
different levels of NAO basis sets (from tier 1 to tier 4) . The reference curve (denoted 
as "ref." ) is obtained with NWChem and a gaussian cc-pV6Z basis, (a): results without 
BSSE corrections; (b): BSSE corrected results. The insets magnify the equilibrium 
region. 



corrected binding energies for tier 2 are in fact almost the same as for tier 4. The latter 
agrees with results from a GTO cc-pV6Z basis (NWChem) within 1-2 meV (almost 
indistinguishable in figure E]). 

Figure [9] demonstrates the same behaviour for a different test case, the binding 
energy of the water dimer using HF (left) and the PBEO [30117] hybrid functional (right). 
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Figure 9. (color online) Convergence of the HF and PBEO binding energies of the 
water dimer (at the PBE geometry, pictured in the inset) as a function of NAO basis 
size (tier 1, 2, 3 for the first three points, and tier 3 for H plus tier 4 for C for the last 
point). Results both with and without counterpoise BSSE correction are shown. The 
dotted line marks the NWChem/cc-pV6Z value. 



The geometry of the water dimer has been optimized with the PBE functional and a 
tier 2 basis. For convenience, the binding energy is computed with reference to H 2 
fragments with fixed geometry as in the dimer, not to fully relaxed H2O monomers. 
This is sufficient for the purpose of the basis convergence test here. Detailed geometrical 
information for water dimer can be found in [01]. In figure [9] the dotted line again marks 
the NWChem GTO cc-pV6Z reference results. Similar to the case of N 2 , we observe 
that the HF binding energy is fairly well converged at the tier 2 level, particularly after 
a counterpoise correction, which gives a binding energy that agrees with the NWChem 
reference value to within 1-2 meV. The BSSE arising from insufficient core description is 
reduced for the PBEO hybrid functional, where only a fraction (1/4) of exact-exchange 
is included. 

In practice, HF calculations at the tier 2 level of our NAO basis sets yield accurate 
results for light elements. Counterpoise corrections help to cancel residual total-energy 
errors arising from a non-HF minimal basis. However, even without such a correction, 
the convergence level is already pretty satisfying (the deviation between the black and 
the red curve in figure [9] is well below 10 meV for tier 2 or higher). 

5.1.2. Heavy elements For heavier elements (Z> 18), the impact of non-HF core basis 
functions on HF total energies is larger. However, the error again largely cancels in 
energy differences, as will be shown below. In order to avoid any secondary effects from 
different scalar-relativistic approximations to the kinetic energy operator, in figure UM 
(upper panels) we first compare the convergence of non-relativistic (NREL) HF total 
energies with NAO basis size for the coinage metal dimers CU2, Ag 2 , and Au 2 at fixed 
binding distance d=2.5 A. (The experimental binding distances are 2.22 A pH], 2.53 A 
[1651 M and 2.47 A [161], respectively.) A gain, we find that KLI-derived minimal basis 
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Figure 10. (color online) Convergence with basis size of the non-relativistic (NREL) 
HF total energies (upper panels), non-relativistic binding energies (middle panels) 
and scalar-relativistic (scaled ZORA 158, 64 ) binding energies (bottom panels) of 
Cu2 (left), Ag 2 (middle), and Au2 (right), at fixed bond length, d=2.5 A. Similar to 
figure results are shown for two sets of minimal basis generated using LDA and 
KLI atomic solvers. For clarity the NREL total energies are offset by -89197.12 eV, 
-282872.42 eV, and -972283.08 eV respectively for Cu 2 , Ag 2 , Au 2 , which correspond 
to the actual values of the last data points with KLI minimal basis. All binding 
energies are BSSE-corrected. The dashed horizontal lines in the bottom panels mark 
the NWChem reference values using aug-cc-pV5Z-PP basis with ECP. 

sets are noticeably better converged (lower total energies) than LDA-derived minimal 
basis sets. In contrast to N2, however, absolute convergence of the total energy is here 
achieved in none of these cases, and the discrepancy increases from Cu (nuclear charge 
Z= 29) to Au (Z= 79). 

For comparison, the middle and lower panels of figure [10] show non-relativistic and 
scalar-relativistic binding energies for all three dimers. The scalar-relativistic treatment 
employed is the scaled ZORA due to Baerends and coworkers |158] (for details of our 
own implementation, see [64]). In all three cases, the binding energies are converged to 
a scale of ~0.02 eV, at least two orders of magnitude better than total energies. In other 
words, any residual convergence error due to the choice of minimal basis (LDA or KLI 
instead of HF) cancel out almost exactly. To compare our prescription to that generally 
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used in the quantum chemstry community where effective core potentials (ECP) are used 
to describe the core electrons and the relativistic effect, we also marked in figure [TU1 the 
reference values computed using NWChem and the aug-cc-pV5Z-PP basis |167[ 1168] . 
The agreement between our all-electron approach and the GTO-ECP one is pretty 
decent for Ag 2 and Au 2 , whereas a larger discrepancy of ~0.03 eV is observed for Cu 2 . 
For this particular case we suspect the remaining disagreement is an issue with respect 
to the atomic reference energy for the Cu atom between both codes. More work will be 
done to fully unravel the point. 

5.2. MP2, RPA, and GW calculations 

In the implementation described here, MP2, RPA, and GW methods require the explicit 
inclusion of unoccupied single-particle states. As a consequence, noticeably larger basis 
sets are needed to obtain converged results in these calculations [169[ I170[ I171[ I172[ 
I173[ 11741 1175[ |38j EE1 EQ]. Much experience has been gained in the quantum chemistry 
community to construct Gaussian basis sets for correlated calculations |157[ 1176] . but 
for NAOs this is not case. In this section, we show how our standard NAO basis sets 
perform for MP2, RPA, and GW calculations, for both light and heavy elements. For 
clarity, we separate the discussions for the convergence of binding energies (in the case of 
MP2 and RPA) and quasiparticle excitations (in the case of GW and MP2 self-energy 
calculations). In contrast to the cases of HF and hybrid density functionals, BSSE 
corrections for RPA and/or MP2 are essential to obtain reliable binding energies. This 
results directly from the larger basis sets required to converge the MP2 or RPA total 
energy pEES HTQl EH EEl EE1 EEJ M\ , yielding larger BSSE for finite basis set size. With 
our standard NAO basis sets, the actual BSSEs in MP2 and RPA (based on the HF 
reference, denoted as RPAOHF in the following) calculations are plotted in figure [11] 
for the example of N 2 . The size of BSSEs in these cases is huge and does not diminish 
even for the pretty large tier 4 basis. It is thus mandatory to correct these errors in 
MP2 and RPA calculations to get reliable binding energies. As one primary interest in 
this work is the applicability of standard NAO basis sets for MP2 and RPA, all binding 
energies presented are therefore counterpoise-corrected. In all HF reference calculations 
in this session, the KLI minimal basis is used. For RPA, GW, and MP2 self-energy 
calculations, we use 40 imaginary frequency points on a modified Gauss- Legendre grid 
(J7j), which ensures a high accuracy for the systems studied here. 

5.2.1. Binding energies As illustrating examples for light elements, in figure IT2l the 
BSSE-corrected MP2 and RPA binding energies for N 2 and the water dimer are shown 
as a function of the NAO basis set size. The dotted line marks reference results computed 
with FHI-aims and the Dunning aug-cc-pV6Z basis. In the case of MP2, the FHI-aims 
aug-cc-pV6Z values agree with that of NWChem to within 0.1 meV. Unfortunately, a 
similar independent reference is not available for RPA, but excellent agreement is also 
seen with smaller basis sets, for which reference RPA data are available for N 2 [38J. 
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Figure 11. (color online) BSSEs in MP2 and RPA@HF binding-energy calculations 
for N2 (d = 1.1 A) as a function of the NAO basis set size. The four points corresponds 
to NAO tier 1 to tier 4 basis sets, respectively. 
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Figure 12. (color online) Convergence of BSSE-corrected MP2 and RPA@HF binding 
energies for N2 (d = 1.1 A) and the water dimer (PBE geometry) as a function of the 
NAO basis set size. The first four points corresponds to NAO tier 1 to tier 4 basis sets 
and the last point corresponds to the composite "tier 4 + a5Z-d" basis. The dotted 
horizontal line marks the aug-cc-pV6Z results. 



Upon increasing the basis size, the biggest improvement occurs when going from tier 
1 to tier 2, with further, smaller improvements from tier 2 to tier 4. For the strongly 
bonded N2 the MP2 binding energy at tier 4 level deviates from the aug-cc-pV6Z result 
by ~ 120 meV, or ~ 1% of the total binding energy. For the hydrogen-bonded water 
dimer, the corresponding values are ~ 3 meV and ~ 1.5% respectively. The convergence 
quality of RPA results with respect to the NAO basis set size is similar. 

Going beyond our FHI-aims standard NAO basis sets, further improvements arise 
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Figure 13. (color online) Convergence of the BSSE-corrected MP2 and RPA@HF 
binding energies for N2 (d = 1.1 A) and the water dimer (PBE geometry) as a function 
of the NAO basis set size. The four points correspond to the tier 1 to tier 4 basis. 
The NAO basis functions are generated for three onset radii (3 A, 4 A, and 6 A) for 
the confining potential. The dotted horizontal line marks the aug-cc-pV6Z results. 

by adding (ad hoc, as a test only) the diffuse functions from a GTO aug-cc-pV5Z basis 
set, denoted "a5Z-d" in the following. The results computed using this composite "tier 
4 + a5Z-d" basis are shown by the last point in figure [12j The deviation between the tier 
4 and aug-cc-pV6Z results is then reduced by more than a factor of two. For the water 
dimer, for example, "tier 4 + a5Z-d" gives -220.9 meV and -206.5 meV for the MP2 
and RPA@HF binding energies, comparable to the quality of the cc-pV6Z basis which 
yields -221.1 meV and -206.9 meV, respectively. Both then agree with the aug-cc-pV6Z 
results (-222.3 meV for MP2 and -208.9 meV for RPA@HF) to within ~ 2 meV. 

In this context, it is interesting to check if the cut-off radii of our NAO functions 
have any influence on the convergence behaviour demonstrated above. As described in 
[64] . the NAO basis functions are strictly localized in a finite spatial area around the 
nuclei, and the extent of this area is controlled by a confining potential. For the default 
settings used in the above calculations, this potential sets in at a distance of 4 A from 
the nucleus and reaches infinity at 6 A. The question is what would happen if we reduce 
or increase the onset radii of this confining potential? The answer to this question is 
illustrated in figure [13] where basis convergence behaviour for N2 and the water dimer 
are shown for three different onset distances of the confining potential. From figure [TBI 
one can see that increasing the onset radius of the confining potential (i.e., enlarging 
the extent of the NAO basis functions) from the default value (4 A) has little effect on 
the convergence behaviour for N 2 or (H 2 0)2- Upon reducing it, noticeable changes of 
the results only occur for tier 1 or tier 2 in certain cases, but the overall effect is very 
small and does not change the general convergence behaviour described above. This 
finding holds in general for covalent and hydrogen bonds. In practice, the onset radius 
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Figure 14. (color online) Convergence of the BSSE-corrected MP2 binding energy 
curve for Au2 with respect to the optimized NAO basis set size (tier 1 to tier 4). 
Results from an independent, Gaussian-type calculation (aug-cc-pV5Z-PP basis set 
with ECP, NWChem code) are included for comparison. 

may always be invoked as an explicit convergence parameter — for instance, much more 
weakly bonded (dispersion bonded) systems benefit from slightly larger radii (5 A - 6 
A) in our experience. Further details on this can be found in |177j . 

We next illustrate the NAO basis convergence for heavy elements, using Au2 as 
an example. In figure HH the MP2 binding energy for the Au2 dimer as a function 
of the bond length is plotted for different NAO basis sets. Relativity is again treated 
at the scaled ZORA level |158[ 164] . The binding curves shown here demonstrate that 
the same qualitative convergence behaviour as for our light-element test cases carries 
over. In essence, significant improvements are gained from tier 1 to tier 2, and basis 
sets between tier 2 and tier 4 yield essentially converged results. For comparison, we 
show a completely independent (NWChem calculations) curve with Gaussian "aug-cc- 
pV5Z-PP" basis sets |f 67] 1168] . The resulting binding energy curve yields rather close 
agreement with our all-electron, NAO basis set results. 

5.2.2. Quasiparticle energies After the discussion of binding energies, we next examine 
how the GW and MP2 quasiparticle energy levels converge with our NAO basis sets. 
The 67 iy°@HF and MP2 quasiparticle HOMO levels for N 2 and the water dimer are 
plotted in figure [15] as a function of basis set size. MP2 is denoted here as "MP2- 
QP" to emphasize that the MP2 self-energy ( |32|) is used, rather than a MP2 total- 
energy difference. We again take the results of an "aug-cc-pV6Z" GTO calculation as 
a reference. The first four data points in each sub-plot correspond to the NAO basis 
sets. The last point represents the composite lu tier 4 + a5Z-d" basis as described 
above. Once again, the biggest improvement occurs when going from tier 1 to tier 2. 
However, the deviation of the HOMO levels from the reference values at tier 2 level is 



40 




Figure 15. (color online) Convergence of the quasiparticle HOMO level of N2 and 
the water dimer obtained with the MP2-QP and G°W°@HF self-energies versus basis 
size. The last data point corresponds to the composite "tier 4 + a5Z-d" basis. The 
aug-cc-pV6Z result is marked by the dotted horizontal line. 



still considerable, ~ 0.2-0.3 eV for G°W° and ~ 0.1-0.2 eV for MP2-QP. These errors 
are brought down to ~0.1 eV for G°W° and ~0.06 eV for MP2-QP by going to a pure 
tier 4 NAO basis set. The remaining error is then further reduced by a factor of two by 
including the diffuse "a5Z-d" part of a 5Z GTO basis set. Accounting for the possible 
underconvergence of the aug-cc-pV6Z itself (compared to the CBS limit), we expect 
an overall ~0.1 eV under- convergence of the composite Ui tier 4 + a5Z-d" results given 
here. This is still an acceptable accuracy, considering the generally known challenge of 
converging the correlation contribution involving virtual states using local orbital basis 
functions PSJ UM EH 0221 ESI EEl ED] ■ Therefore 111 tier 4 + a5Z-d" basis sets were 
used in the benchmark calculations presented in the section |6j 

5.2.3. G°W° calculations for the benzene molecule G°W° calculations for molecules 
have been reported in a number of publications in recent years using various numerical 
frameworks and computer code packages. In particular in the solid-state community, 
calculations have been based on plane wave methods together with the supercell 
approach and pseudopotential approximations. Their advantage is a systematically 
convergable basis set (plane waves) especially for the unoccupied spectrum. However, 
the two approximations mentioned above (pseudopotential and supercell) can be drastic 
[1751 11751 11501 H5T1 IT821 H551 H53] . In particular, the Coulomb operator in vacuum is 
not screened in a standard plane wave approach. As a result, different images of the 
system interact with one another across supercell boundaries |180[ I181[ I182[ I183[ 1184] . 
In addition, the slow converence of G°W° results with the plane-wave cutoff of the 
unoccupied spectrum is notorious |175[ I185[ 1186] . 

Table [3] reports literature G°W° results for the benzene molecule |187[ I128[ I188[ 
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Table 3. G°W° HOMO and LUMO values for the benzene (C 6 H 6 ) molecule 
obtained by FHI-aims and several other numerical approaches as reported in literature. 
"A.E" and "P.P." in the second column refer to "all electron" and "pseudopotential" 
respectively. 



G°lT -type 


P.P./A.E. 


basis type 


HOMO (eV) 


LUMO (eV) 






NAO "tier 3" a 


-8.99 


1.06 




A.E. 


AT A /~\ U ± " A )) 17, 

NAO tier 4 


-9.06 


0.96 






NAO "tier 4 + a5Z-d" a 


-9.05 


0.94 


G°W°@LBA 


P.P. 


plane waves + NAOs b 
plane waves + extrapol. c 


-9.03 
-9.10 


1.54 

/ 




plane waves + Lancos d 


-9.40 


/ 






NAO "TZDP" e 


-8.78 


1.24 






real-space grid-^ 


-9.88 


0.47 


G°W°@RF 


A.E. 


NAO "tier 4" a 
Gaussian "cc-pVTZ" 9 


-9.64 
-9.28 


1.51 

/ 


Experiments 






-9.24 h 


1.12 1 



"this work 
b Ref. [TEH] 
c Ref. [HI] 
d Rcf. [188] 
e Ref. H53| 
,l Ref. USD] : 

*Re£ W ■ 



plane wave basis augmented with SiESTA-type localized atomic orbitals 
plane wave basis extrapolated to infinite energy cutoff 
plane wave basis plus Lanczos trick to remove the empty states 
'Ref. OH] 9Ref. [US] 

(negative) ionization energy (IE). The vertical IE only differs from 
this value by 0.01 eV according to the NIST database [191J . 
(negative) vertical electron affinity (IE) 



186, I153[ I189j . a particularly often studied case, in comparison to our own results. 
We focus on G°W°@LDA and G°W°@RF for the HOMO and LUMO levels. Based 
on these results, it is clear that there is a significant degree of scatter, even between 
results that are ostensibly converged using the same fundamental approximations. Our 
own results for NAO "tier 3" , "tier 4" , and "tier 4 + a5Z-d" basis sets suggest internal 
convergence at the "tier 4" level: -9.06 eV for the HOMO, and 0.94 eV for the LUMO in 
G°jy°@LDA, compared to -9.64 eV and 1.51 eV in G W°@RF. The LUMO values are 
unbound and can be interpreted as experimental resonance energies (here taken from 
the tabulated negative vertical electron affinity). In either case (HOMO or LUMO, 
G°W°@LDA or G W°©HF) the results are not far from the experimental values. 

The same cannot be said for the comparison between the different numerical 
implementations. For instance, the HOMO values in G°W°@LDA range from —8.78 eV 
(small numeric basis set and pseudopotentials) to —9.88 eV on a real-space grid. Even 
within the plane wave based approaches, the values range from —9.03 eV to —9.40 eV. 
Similar discrepancies arise for G°W°@BF, and for the LUMO values. Within the scatter 
evidenced by Tabled we believe that our own results possess some merit, as the system 
(i) is isolated (no supercell), (ii) is treated fully all-electron, and (iii) any residual 
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convergence issue of the unoccupied spectrum should be exposed by the diffuse Gaussian 
basis functions ("a5Z-d"). 

6. Basis set converged benchmark data for the G2 and S22 molecular test 
sets 

Our final section utilizes the preceding methodologies and techniques to provide accurate 
all-electron results for molecular test sets close to basis set convergence. We cover vertical 
ionization energies and binding energies at different levels of theory, for well-defined, 
published molecular geometries. 

6.1. Vertical ionization energies from HF, MP 2, and G°W° methods 

The G°W° method has been used to calculate the single-particle properties for various 
small- and medium-size molecules with considerable success |1931 11941 11951 H87[ I196j 11971 
IT981 IT991 1200] lEEJ EDB EQ2J QS21 [153]. The dependence of G°W° on the starting point 
has also been looked at in the past [601 ED H97[ [202J. Here, using a selected collection 
of atoms and molecules from the G2 ion test set for ionization energies [79], we aim 
to establish the overall performance of G°W° for computing IEs close to the basis set 
limit, and systematically examine its dependence on the starting point. The selection 
of molecules is based on the availability of experimental geometries and experimental 
vertial IEs. With vertical IEs we denote ionization energies at fixed geomerty, i.e. no 
structural relaxation after excitation. This quantity is directly comparable to the GW 
and MP2-QP and quasiparticle energies. We will also assess the MP2-QP approach for 
determining IEs since this method was used in quantum chemistry in past decades, but 
direct comparisons in the literature are scarce [H] . Finally, results for IEs determined by 
MP2 total energy differences (denoted here simply as "MP2" in contrast to "MP2-QP" ) 
between the neutral atoms/molecules and the corresponding positively charged ions are 
also presented. We note that IEs given by MP2-QP essentially correspond to the MP2 
total energy difference if the HF orbitals of the neutral systems were used in the ionic 
calculation (for a discussion, see, e.g., [B]). All calculations were carried out at the 
experimental geometries. As mentioned above, the composite ll tier 4 + a5Z-d" basis 
set is used for most of the elements except for a few cases (e.g., rare gas atoms) where 
the tier 4 NAO basis is not available. In these cases the full aug-cc-pV5Z functions 
are added to the NAO tier 2 basis. On average we expect the chosen basis setup to 
guarantee a convergence of the HOMO quasiparticle levels to ~0.1 eV or better for the 
calculations presented in the following. 

In figure [TBI we present histograms of the error distributions given by HF, MP2-QP, 
G°W°@B.F and G°W°@PBE0 for this database. The actual IE values are presented in 
[7J Compared to HF, one can see that the deviations from the experimental values are 
much smaller in G°W° and MP2-QP. On average HF tends to overestimate IEs. This 
trend is corrected by MP2 and G°W°. Using the same HF reference, the magnitude of 
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Table 4. Mean error (ME), mean absolute error (MAE), and mean absolute percentage 
error (MAPE) for the ionization energies in the G2-I subset computed with HF, MP2- 
QP, MP2, and G°W° based on HF, PBE, and PBEO. 
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Figure 16. Histograms of the error distribution for IEs for a set of 50 atoms and 
molecules calculated with HF, G°W°@HF, G°W°@PBE, G°VF°@PBE0, MP2-QP, and 
MP 2. 

the correction is smaller in G°W° than in MP2-QP, due to the renormalization effect 
coming from the screened Coulomb interaction in the GW self-energy. Concerning the 
starting-point dependence of the G°W° method, G°W° based on HF gives too large IEs 
on average, whereas G°W° based on PBE does the opposite. Consequently, G°W° based 
on the hybrid density functional PBEO appears to be the best compromise, although 
a slight underestimation of the IEs is still visible. Furthermore, comparing MP2-QP 
with MP2 shows that the two approaches yields very similar results, implying that 
for the light elements reported here orbital relaxation effects are not significant [203J. 
Table H] summarizes the error statistics. For the subset of atoms and molecules MP2-QP 
gives the smallest mean error (ME) and G°W°@PBE0 the smallest mean absolute error 
(MAE). 
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6.2. Benchmark MP2 and RPA results for the G2-I atomization energies 

The 55 atomization energies from the original G2-I set (Table III of |204j ) serve as a 
good benchmark database for ground-state total-energy methods for covalently bonded 
systems. In Table |5] we present our MP2 and RPA atomization energies for the G2-I set 
as computed using FHI-aims and the "tier 4 + a5Z-d" basis set. The geometries used 
in the calculations are those determined by Curtiss et al. [204\, i.e., all-electron MP2 
optimization with a Gaussian 6-31G* basis set. The reference data, taken from [205J, is 
derived from experiment and corrected for zero-point vibrations. 

Our MP2 and RPA results are in good agreement with those reported by Feller 
and Peterson [206J and Paier et al. [90], respectively, both obtained with GTO basis 
sets and extrapolated to the complete basis set limit. Table [5] demonstrates that MP2 
performs best for these small covalently bonded molecules. The RPA approach, on 
the other hand, which has a broader applicability (e.g. bond breaking and/or metals 
where MP2 fails), exhibits significant underbinding. For RPA we also observe that KS- 
PBE provides a much better reference than HF. Recently computational approaches to 
overcome the underbinding behaviour of the standard RPA scheme have been proposed. 
These comprise contributions from single excitations [ID] and second-order screened 
exchange pTOl [90]. 

6.3. Benchmark MP2 and RPA binding energies for the S22 molecular set 

The S22 molecular set proposed by Jurecka et al. [80J has become a standard benchmark 
database for testing the accuracy of existing and newly developed methods for the 
description of weak interactions. S22 represents an "unbiased" set in the sense that it 
contains molecules of different bonding nature (7 hydrogen bonded, 8 dispersion bonded, 
and 7 with mixed bonding character) and of different size (ranging from small ones like 
the water dimer to relatively large ones like the Adenine-thymine dimer containing 30 
atoms). The MP2 and CCSD(T) interaction energies for this set of molecules were 
already computed by Jurecka et al. and extrapolated to the Gaussian complete basis 
set (CBS) limit. A more consistent and accurate extrapolation for the CCSD(T) values 
was recently carried out by Takatani et al. [207] . which we therefore adopt as reference 
here. While MP2 calculations for S22 are common, RPA benchmark calculations for the 
whole set have not been performed, yet. RPA-type calculations for S22 have recently 
been reported [98 1 l4"0~ t [208] . The agreement between the data from different groups is not 
perfect, and the basis incompleteness could be an issue. Of our own RPA numbers only 
the MAEs and MAPEs have been reported previously [10] . Now in Table [6] the actual 
MP2, RPA@HF, and RPA@PBE binding energies for these molecules as computed using 
FHI-aims and the composite "tier 4 + a5Z-d" basis set. With our basis setup, the MP2 
binding energies are underestimated by ~4.5 meV (2% on a relative scale) compared 
to the MP2/CBS results reported in [80]. We expect a similar convergence of the 
RPA numbers based on our basis set convergence tests shown in previous sections. 
Compared to the CCSD(T) reference data, RPA@PBE, which nowadays dominates 
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Table 5. G2-I atomization energies computed with MP2, RPAQHF, and RPAQPBE. 
All numbers have been calculated with FHI-aims and "tier 4 + a5Z-d" basis sets. The 
experimental reference results are taken from [205 . 



A/Tolocn 1p 


MP2 


RPA9HF 

J.V.A -Tl. :\ 111 


RPAOPRE 




BeH 


54.3 


55.6 


49.9 


49.8 


C2H2 


412.9 


374.1 


377.3 


403.4 


C2H4 


564.3 


527.4 


533.1 


562.4 


C2H6 


709.9 


671.4 


678.7 


712.7 


CH 


80.0 


75.7 


80.5 


83.9 


CH2( 1 A 1 ) 


176.6 


168.1 


173.3 


180.9 


CH2( 3 B 3 ) 


189.3 


181.0 


178.2 


190.2 


CH 3 


304.6 


290.5 


292.4 


307.9 


CH3CI 


397.1 


370.1 


368.0 


394.5 


CH3OH 


515.7 


475.8 


486.2 


512.8 


CH 3 SH 


466.5 


432.2 


446.2 


473.7 


CH 4 


416.1 


395.8 


402.4 


420.1 


CN 


165.8 


138.6 


170.1 


181.0 


CO 


270.8 


234.7 


241.9 


259.3 


co 2 


412.1 


350.9 


358.9 


389.1 


cs 


178.3 


147.9 


154.7 


171.2 


Cl 2 


64.5 


50.1 


45.8 


58.0 


C1F 


67.2 


47.2 


47.4 


61.5 


CIO 


61.3 


43.3 


56.7 


64.6 


F 2 


42.6 


16.0 


28.4 


38.2 


H2CO 


381.4 


343.1 


351.6 


373.7 


H20 


235.6 


212.1 


220.8 


232.7 


H202 


274.4 


231.7 


252.2 


268.7 


H2S 


179.8 


167.6 


170.0 


182.5 


HCN 


322.3 


281.4 


295.6 


311.7 


HCO 


287.2 


251.3 


260.3 


278.5 


HC1 


107.0 


98.0 


92.7 


106.4 


HF 


145.2 


129.4 


130.3 


141.1 


HOC1 


168.5 


138.3 


148.8 


164.6 


Li 2 


18.2 


13.5 


18.3 


24.4 


LiF 


145.2 


127.3 


125.4 


138.9 


LiH 


53.2 


50.5 


54.1 


58.0 


->2 


OQQ n 

zoo.U 


lyo.o 


01 n 
ziy.y 


zzo.o 


N 2 H 4 


436.4 


393.8 


421.5 


438.6 


NH 


79.2 


74.1 


81.3 


83.6 


NH 2 


178.0 


164.6 


177.3 


182.0 


NH 3 


295.3 


271.9 


288.0 


298.0 


NO 


155.2 


120.2 


145.5 


152.8 


Na 2 


13.1 


8.7 


14.2 


17.0 


NaCl 


101.3 


91.5 


88.7 


97.8 


02 


129.8 


95.4 


110.8 


120.3 


OH 


106.6 


96.5 


102.0 


106.7 


P 2 


118.1 


92.2 


112.1 


117.2 
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Table \E[ (continued) 


Molecule 


MP 2 RPA@HF 




i _ . . 

hjxp. 


PH 2 


145. 1 


140. / 


148.2 


153.1 


PH 3 


oo i a 

231.0 


221. / 


231.9 


243.6 


s 2 


100.6 


79.1 


89.9 


101.7 


SO 


132.7 


104.6 


110.9 


125.0 


so 2 


273.8 


206.6 


232.9 


258.4 


c ■ 

012 


73.7 


64.6 


67.7 


74.7 




518.4 


500.9 


506.6 


530.6 


SiH2( 1 A 1 ) 


146.4 


140.8 


146.5 


151.7 


SiH2( 3 Si) 


127.0 


124.0 


124.7 


130.9 


SiH 3 


219.2 


213.4 


217.1 


227.0 


SiH 4 


314.5 


306.2 


310.4 


322.0 


SiO 


204.3 


170.0 


179.0 


191.7 


ME 


1.0 


-21.5 


-13.3 




MAE 


5.9 


21.7 


13.3 




MAPE 


4.3 % 


13.4% 


7.6% 





RPA-type production calculations, systematically underestimates all of the three weak 
bonding categories and gives a MAE of 39 meV and a MAPE of 16%. If instead HF is 
used as a starting point (i.e., as for MP2), the description of hydrogen bonding improves, 
while the description of dispersion bonding worsens. The overall MAE for RPA@HF is 
51 meV and the MAPE 25%. We are thus faced with the conundrum that RPA can 
describe the weak interactions that are beyond the reach of LDA, GGA, and hybrid 
functionals, but the accuracy of the two standard RPA schemes is not spectacular. We 
have analyzed the origin of this underbinding behaviour in [40] and proposed a simple 
solution to overcome this problem. 

7. Conclusions and outlook 

To summarize, we have presented a resolution of identity framework for the two-electron 
Coulomb operator that allows efficient, accurate electronic structure computations based 
on HF, hybrid functionals, MP2, RPA, and GW based on the flexible basis function 
form of NAOs. We have shown that NAO basis sets as implemented in FHI-aims 
are a competitive choice for approaches involving exact-exchange and/or non-local 
correlation terms, with rather compact basis sets sufficing for essentially converged 
results. Our simple "on-the-fly" scheme to construct the auxiliary basis functions using 
Gram-Schmidt orthonormalization of the "on-site" products of single-electron atomic 
orbitals gives a natural, accurate representation of the two-electron Coulomb operator 
for practical calculations. Taken together, our framework paves the way for an extended 
usage of NAOs in more advanced computational approaches that go beyond LDA and 
GGAs. Specifically, we have applied the G°W° and MP2 quasiparticle approaches to 
compute the vertical IEs of a set of small molecules, and the RPA method to compute 
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Table 6. Binding energies (in meV) of the S22 molecular set [50] calculated with MP2, 
RPA@HF, and CCSD(T). MP2 and RPA results have been obtained with FHI-aims 
and "iier 4 + a5Z-d" basis sets. The reference CCSD(T) results are from [207] . 





Molecules 


MP 2 


RPA@HF 


RPAQPBE 


CCSD(T) 


1 


NH3 dimer 


136 


116 


112 


137 


2 


H2O dimer 


214 


201 


182 


218 


3 


Formic acid dimer 


799 


790 


744 


816 


4 


Formamide dimer 


681 


662 


645 


700 


5 


Uracil dimer 


881 


846 


816 


898 


6 


2-pyridoxine . 2-pyridoxine 


749 


666 


678 


738 


7 


Adenine thymine 


714 


639 


658 


727 


8 


CH4 dimer 


21 


11 


17 


23 


9 


C2H4 dimer 


67 


37 


49 


65 


10 


Benzene. CH4 


78 


35 


49 


63 


11 


Benzene dimer 


211 


40 


82 


114 




(slip parallel) 










12 


Pyrazine dimer 


296 


110 


144 


182 




Uracil dimer 


ARO 


olU 


/ y 




1 1 


Indole. benzene 


040 


^f; 

oO 


1 A R 


1 QQ 
lyy 


15 


Adenine. thymine (stack) 


ti A 1 

641 


388 


A OO 

422 


506 


l(i 


Ethene.ethine 


72 


57 


54 


65 


17 


Benzene. H2O 


153 


119 


123 


143 


18 


Benzene. NH3 


114 


74 


85 


101 


19 


Benzene. HCN 


223 


178 


170 


197 


20 


Benzene dimer (T-shape) 


156 


75 


96 


118 


21 


Indole. benzene (T-shape) 


200 


182 


215 


244 


22 


Phenol dimer 


334 


250 


266 


308 




ME 


26 


51 


39 






MAE 


37 


51 


39 






MAPE 


19% 


25% 


16% 





the G2-I atomization energies and the interaction energies for the S22 molecular set. We 
believe that the well converged numbers reported in this work may serve as benchmarks 
for future studies. Beyond the specific examples given here, our RI framework as a 
whole has already proven to be stable and mature in a number of scientific applications 

[9mi4^i4ii [iTn[Mi2iBim] . 

Based on the results above, NAOs emerge as a promising route towards more 
compact basis sets for correlated methods. Among our ongoing developments is 
the attempt to design more optimized NAO basis sets for MP2, RPA, and GW 
calculations. Another active line of development is to improve the scaling behaviour 
of the computational cost with system size, exploiting the locality of the NAO basis 
functions |147] . Last but not least, we are working on the extension of the present 
scheme to periodic systems. All this combined, we expect NAO-based implementations 
of methodologies that go beyond LDA and GGAs to become competitive alternatives 
to the traditional implementations that are based on GTO or plane-wave basis sets, in 
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particular towards the limit of very large systems: the benefit of compact basis set size 
at a given accuracy level should help tackle problem sizes that, otherwise, could not be 
done. 
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Appendix A. Matrix elements for numeric atom-centered orbitals 

Appendix A.l. Coulomb potential of a numerical radial function 

To reduce the (formally) six-dimensional Coulomb integrations to three dimensional 
ones, we first solve Poisson's equation for each P M (r') (classical electrostatics). We 
define Q^ijc) as 

Q» = J drXr-r')P^(r'). (.1) 
Using the Laplace expansion of the Coulomb potential v (r — r') = l/|r — r'|, 

A-jr if ^ 

v ( r - r ') = E ^l Y im(m*m(r') x pfj, (.2) 

lm < 

where r< = min(r, r') and r > = max(r, r'), the integration of the angular part in (TTJ) 
can be done analytically, yielding 

Q„(r) = ^l Ylm (h). (.3) 

The radial part a a i K (r) is given by a simple one-dimensional (numerical) integration 
[2T21 1551 IM] 



An 



21 + 1 



(-4) 



Thus, the three-center and two-center Coulomb integrals fl37j) and ff45|) reduce to the 
three-dimensional integrals in ( 1661) and ( 1671) . 

Appendix A. 2. Grid-based three-center and two-center integrals 

The three-center and two-center integrals in the present work are computed by grid- 
based integrations using overlapping, atom-centered spherical grids and the same 
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technology that is used in many quantum-chemical applications for the exchange- 
correlation matrix in DFT |154[ 165] . The integration grid points r = r(a, s, t) are 
uniquely specified by the atomic center a, the radial shell number s, and angular point 
t. w(s,t) is the corresponding integration weight. Details of our own implementation 
(FHI-aims) are given in [64"1 1155] . Since these are true three-center integrals, we restrict 
the integration domain for a particular integral element to the grids associated with the 
atoms on which the basis functions in question are centered. For instance, denoting the 
respective atoms by a±, a 2 , 03, the three-center integrals in ( 1661) can then be discretized 
as 

(ij = Yl ^^(a, r)«?(s,t)0i(r)0j(r)Q M (r), (.5) 

where ^3(0, r) is a three-center partition function that satisfies 

Y Ps(a,r) = l 

a=ai ,a,2> a 3 

everywhere in the overlapping region of the three functions, and is zero otherwise. The 
underlying numerical grids can in principle be increased up to arbitrary accuracy if 
needed. The two-center integrals ( 1671) can be performed in a similar fashion using 
overlapping grids, or with the spherical Bessel transform techniques explained below. 

Appendix A. 3. Two-center integration in Fourier space 

As mentioned in section 14.31 two-center integrals of numeric atom-centered basis 
functions like in (1371) and S M „ in ( 1411 can be efficiently calculated in Fourier space 
as described by Talman [1321 1213] . This and the following flTJ subsections give the 
details of our implementation. We first describe the general procedure here, and the 
central ingredient of our implementation, the logarithmic spherical Bessel transform 
(logarithmic SBT, logSBT), will be presented in the next subsection. 

As is well-known, the overlap of two functions f(r) and g(r) can be expressed in 
Fourier space as 

J f(r)g(r - R)dr = J f{k)g{-k)e^ R dk. (.6) 

The Fourier transform f(k) in (TjS]) of an atomic function f(r) = f (r)Yi m (f) has the 
same angular momentum in real and Fourier space for symmetry reasons 

f(k) = (27T)-* J f(r)e-* k r dr = r l f~(k)Y lm (k). (.7) 
The radial part f(k) is given by the SBT of f(r) 

Rk) = \Jl 3i{kr)f{r)r 2 dr. (.8) 

If f{r) is tabulated on a logarithmic grid, its SBT can be calculated efficiently on an 
equivalent logarithmic grid using the fast logSBT algorithm |214[ I215[ 1216] as described 
in the next section. 
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The three-dimensional integral fljB]) can be separated by expanding the plane wave 
e ik R j n gphgricgj harmonics and spherical Bessel functions 



- 00 12 

e M H2^)^ ii V^(^)E^W y ^). (-9) 

L=0 M 

The separation yields 

/ f(r)g(r - R) dr = (27r)f ^2r w ' +L I L {R)A L (R) (.10) 

J L 

with a radial integral 

h(R) = ^1 j L (kR)f(k)g(k) k 2 dk (.11) 

and an angular integral 

A L (R) = Y J Y *LM{R)C{lm ] l'm' ] LM). (.12) 

M 

The triple- F integrals 

C(lm;l'm';LM) := / Y lm {k)Y Vm ,{k)Y LM {k) dk (.13) 

in (I.12p can be calculated efficiently using recursion formulae |213j . They are nonzero 
only for L = \l — — + 2, . . . , (l + l'). If two atom-centered functions do not overlap, 
the overlap integrals Il{R) of course vanish. 

For any given distance R, the radial integrals Il{R) in (l.lip can be calculated 
directly using the trapezoidal rule on the logarithmic grid when j^kR) is evaluated 
in logarithmic Fourier space as described in the next section. If integrals of the same 
atom-centered functions for many differing distances are needed, one can compute these 
more efficiently by interpreting (I. lip as an SBT of P(k) = f(k)g(k). By applying the 
logSBT, and interpolate for all needed distances R, one can obtain all the integrals 
at tight-binding cost, meaning that efficient recursion formula (together with spline 
evaluations) can be used intead of evaluating each integral numerically. 

Coulomb interactions of atomic functions can be calculated with comparable ease 
where the integrand in (I. lip is multiplied with the Coulomb kernel Ait/k 2 

2 r, , LD x/(W) ,2 



V L (R)=A7r^-j 3L (kR) JK ^ ' k 2 dk. (.14) 

The Coulomb interaction generally does not vanish even if the two charge densities do 
not overlap. However, it has a simple multipolar behaviour and explicit integration 
of (1.141) can thus be avoided. The function Vl(R) can formally be interpreted as the far 
field of a charge distribution of angular momentum L whose radial part P(r) is given 



51 



by its SBT P{k) = f(k)g(k). Therefore, it only depends on the multipole moment of 
P, which is encoded in its limiting behaviour for small k. From this, it can be shown 
that Vl{R) vanishes for all L < I + I'. For L = I + I' we get 

t/ (m 4tt [2 (2/ + 21'-!)!! 

with (2n — 1)!! = 1 • 3 • • • (2n — 1) and with multipole moments 

1 f°° 

Pf = zT+i /o r2+ ' /(r)dr - ( - 16) 

Therefore, the Coulomb interaction of non-overlapping functions depends only on the 
product of their multipole moments and can also be calculated at tight-binding cost. 

As a side remark we emphasize that two ABFs do not interact via the Coulomb 
interaction if they do not overlap and at least one of the two multipole moments is zero. 
As shown by Betzinger in |217j . all but one of the ABFs for a given atom and a given 
angular momentum Im can then be chosen to be multipole free by means of a suitable 
unitary transformation. 

Appendix A. 4- Logarithmic spherical Bessel transform 

This section describes our implementation of the logSBT algorithm. For a more 
extensive description we refer the interested reader to Talman [216] and Hamilton |215] . 

The SBT as defined in ([78]) can be written as the integral over a kernel J~(kr) and 
a right-hand-side J-(r): 

f(k) = k~ a I a - ^- 3l {kr){krT r^m- (-17) 

= :J(kr) V ' 







The choice of the power bias parameter < a < 3 is crucial for the numerical accuracy 
and stability of this method and will be discussed further below. The basic idea of the 
logSBT is that in logarithmic coordinates (p = log r and k = log k) the kernel reads 
J{k + p) and ( 1.1 7ft turns into a convolution, which can be efficiently calculated using 
fast Fourier transorms (FFTs). Please note that dr jr = dp in (I.17p . 

As pointed out by Hamilton |215j . this procedure is exact if both ^(p) = 
e^ 3 ~ a ^ p f(e p ) and the corresponding SBT term ^F(k) = e aK f(e K ) are periodic in 
logarithmic space and analytic expressions for the logarithmic Fourier transform of the 
kernel are used. Periodicity can be achieved, e.g., by choosing a near 1.5 and using 
a sufficiently wide logarithmic grid. Under these circumstances, both ^(p) and T(k) 
smoothely drop to zero on both ends, which therefore can safely be connected. 

Unfortunately, the scaling factor k~ a turns out to be quite problematic. By design 
of the algorithm, the absolute error of F(k) = k a f(k) before the final scaling is typically 
of the order of machine precision, i.e., about 10 -15 . This is true even if the magnitude 
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of the exact value is much smaller than that. After scaling, however, the absolute error 
can get arbitrarily large because k can be very small on a wide logarithmic grid. 

Talman |216] circumvents this problem by using two separate a for small and large 
k and joining the two results where they differ least. The small-a (a = 0) calculation 
cannot be done assuming periodicity for I = because J-'(p) does not decay to zero for 
p — > — oo. Therefore, a trapezoidal rule is used for the integration, which works well for 
small k where the Bessel function is smooth on a logarithmic scale. This does not break 
with the spirit of logSBT, because the trapezoidal rule can be formulated using FFTs, 
too. 

In order to avoid the second transform, we take a different approach. In practice, 
one needs the SBT only for one kind of integral and it is sufficient to calculate k a f(k) 
to high absolute accuracy for a single given a, which can be used as power bias for the 
transform. 

This works well for all cases but a = and 1 = 0. Here, we cannot simply resort to 
the trapezoidal rule because it is invalid for high k where the Bessel function oscillates 
rapidly. Instead, we separate jo(e T ) into a smooth part proportional to erfc(r/Ar ) and 
a properly decaying rest. We use the first part for a trapezoidal rule and the second 
part for the log-periodic algorithm. Fortunatly, these two schemes differ only in the way 
the kernel is constructed so that the two kernels can simply be added up to a "hybrid" 
kernel. The actual transform is not affected and numerically not more expensive than 
an ordinary logSBT. 

Just like Talman |216] we double the domain during the transforms for / = and 
I = 1 in order to avoid the need for large domains for a proper decay behaviour. 

Appendix B. Ionization energies of the a set of atoms and molecules 

In Table [1] the individual numbers for the vertical ionization potentials for 50 atoms and 
molecules as computed with 6 different computational approaches are presented. The 
calculations are performed with FHI-aims and "tter+a5Z-d" basis set. 

Appendix C. Modified Gauss-Legendre grid 

For the integrals over the imaginary frequency axis (e.g., for the RPA correlation energy, 
equation ( 1601) ), we use a modified Gauss-Legendre quadrature. The Gauss-Legendre 
quadrature provides a way to numerically evaluate an integral on the interval [—1:1] 



where x« and tOj are the integration points and the corresponding weights, respectively. 
For our purposes a transformation procedure is applied to map the integration range 
from [—1 : 1] to [0 : oo] whereby the and Wi have to be changed accordingly. 




(.18) 
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Table 1. The vertical ionization potentials (in eV) for 50 atoms and molecules 
(taken from G2 ion test set [75]) calculated with HF, MP 2, MP2-QP, G°W°@HF, 
G°T4 7 °@PBE0, and G°1F°@PBE in comparison to the experimental values, taken from 
the the NIST database [191] . The mean absolute errors (MAE) for the three approaches 
are also shown. 

Molecule Exp. HF MP2 MP2-QP G°W°@BF G°W°@PBE0 G°W°@PBE 



Al 


5.98 


5.95 


5.85 


6.08 


6.24 


5.94 


5.64 


Ar 


15.76 


16.08 


15.87 


15.61 


16.08 


15.51 


15.21 


B 


8.30 


8.68 


8.33 


8.79 


8.73 


8.11 


7.65 


BClq 


11.64 


12.48 


12.58 


11.65 


12.37 


11.63 


11.25 


BFq 


15.96 


18.04 


16.19 


15.02 


16.77 


15.79 


15.21 


Be 


9.32 


8.48 


8.87 


8.98 


9.16 


9.26 


9.03 


c 


11.26 


11.95 


11.33 


11.79 


11.68 


10.93 


10.47 


C2H2 


11.49 


11.19 


11.75 


11.43 


11.76 


11.29 


11.01 


C2 H4 


10.68 


10.23 


10.77 


10.37 


10.83 


10.47 


10.22 


C 2 H4 S 


9.05 


9.43 


9.27 


8.74 


9.44 


8.94 


8.71 


C9H«;OH 


10.64 


12.05 


11.30 


10.08 


11.47 


10.63 


10.20 




9.25 


9.15 


9.88 


9.08 


9.63 


9.20 


9.00 


y • L L z. -- - 1. 


10.20 


10.31 


10.55 


9.94 


10.70 


10.12 


9.85 


CH 2 S 


9.38 


8.24 


9.76 


8.04 


8.22 


9.27 


9.01 




9.84 


10.47 


9.78 


10.61 


10.28 


9.59 


9.24 


CH3CI 


11.29 


11.87 


11.63 


11.20 


11.88 


11.31 


11.03 


CH?F 


13.04 


14.46 


14.18 


12.98 


14.15 


13.28 


12.77 


CH3SH 


9.44 


9.67 


9.56 


9.25 


9.81 


9.31 


9.06 


CH 4 


13.60 


14.85 


14.46 


14.18 


14.89 


14.27 


13.98 


CHO 


9.31 


11.13 


9.36 


9.92 


10.64 


9.61 


9.14 


CO 


14.01 


15.10 


14.55 


14.06 


14.85 


13.78 


13.30 


C0 2 


13.78 


14.83 


14.86 


13.29 


14.48 


13.68 


13.21 


CS9 


10.09 


10.14 


10.86 


10.03 


10.55 


10.02 


9.72 


CI 


12.97 


13.09 


12.90 


13.21 


13.32 


12.83 


12.51 


CL 


11.49 


12.08 


11.66 


11.38 


12.13 


11.49 


11.04 


C1F 


12.77 


11.63 


11.34 


11.08 


11.72 


11.15 


10.72 


F 


17.42 


18.50 


17.42 


17.30 


17.73 


17.07 


16.71 


FH 


16.12 


17.70 


16.47 


14.93 


16.39 


15.83 


15.39 


H 


13.61 


13.61 


13.61 


13.61 


13.61 


13.04 


12.52 


He 


24.59 


24.98 


24.41 


24.58 


24.68 


24.01 


23.59 


Li 


5.39 


5.34 


5.38 


5.38 


5.68 


5.84 


5.67 


Mg 


7.65 


6.88 


7.40 


7.44 


7.56 


7.64 


7.71 


N 


14.54 


15.54 


14.66 


14.98 


14.85 


14.06 


13.51 


N 2 


15.58 


16.71 


15.48 


17.22 


17.27 


15.45 


14.86 


NH 3 


10.82 


11.70 


11.06 


10.29 


11.36 


10.70 


10.32 


Na 


5.14 


4.97 


5.11 


5.09 


5.37 


5.51 


5.51 


NaCl 


9.80 


9.68 


9.42 


9.98 


9.59 


9.09 


8.79 


Ne 


21.56 


23.14 


21.63 


20.22 


21.76 


21.10 


20.54 





13.61 


14.20 


13.48 


14.64 


13.90 


13.37 


13.04 


o 2 


12.30 


15.22 


11.85 


12.57 


13.71 


12.33 


11.68 


ocs 


11.19 


11.47 


11.95 


11.16 


11.70 


11.16 


10.88 


OH 


13.02 


13.98 


13.14 


13.21 


13.38 


12.79 


12.41 


P 


10.49 


10.67 


10.56 


10.78 


10.72 


10.24 


9.94 


P 2 


10.62 


10.10 


10.86 


10.59 


10.70 


10.35 


10.13 
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Table ( continued ) 



Molecule 



Exp. 



HF 



MP2 MP2-QP G°W°@BF G°W°@PBE0 G°W°@PBE 



PH 3 

S 

s 2 

SH 2 
Si 



10.59 
10.36 

9.55 
10.50 

8.15 
12.30 



10.58 
10.33 
10.38 
10.49 
8.20 
13.24 



10.58 
10.15 

9.38 
10.53 

8.10 
12.81 



10.47 
11.04 

9.74 
10.34 

8.33 
12.90 



10.90 
10.69 
10.21 
10.74 
8.38 
13.33 



10.44 
10.31 

9.39 
10.27 

8.01 
12.68 



10.22 
10.12 

9.05 
10.06 

7.76 
12.29 



SiH 4 



Specifically, we use the modification proposed for the evaluation of the Casimir-Polder 
integral [21"8] : 

% = x£ +Xi)ni - Xi \ (.19) 
with xq set to 0.5. The weights for the tranformed grid are then given by 



This modified Gauss-Legendre scheme allows a quick convergence of the frequency 
integration with a relatively small number of frequency points. In our implementation, 
a 40-point grid gives micro-Hartree total energy accuracy for the systems investigated 
in this work. 
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